! This file is AUTOGENERATED
! To change, edit m_wcml_molecule.m4 and regenerate

module m_wcml_molecule

  use fox_m_fsys_realtypes, only: sp, dp
  use FoX_wxml, only: xmlf_t

#ifndef DUMMYLIB
  use fox_m_fsys_format, only: str
  use m_common_error, only: FoX_error
  use FoX_wxml, only: xml_NewElement, xml_EndElement
  use FoX_wxml, only: xml_AddAttribute, xml_AddCharacters, xml_AddNewline

! Fix for pgi, requires this explicitly:
  use m_wxml_overloads
#endif

  implicit none
  private

  interface cmlAddMolecule
    module procedure cmlAddMoleculeSP
    module procedure cmlAddMoleculeSP_sh
    module procedure cmlAddMolecule_3_SP
    module procedure cmlAddMolecule_3_SP_sh
    module procedure cmlAddMoleculeDP
    module procedure cmlAddMoleculeDP_sh
    module procedure cmlAddMolecule_3_DP
    module procedure cmlAddMolecule_3_DP_sh
  end interface

  interface cmlAddAtoms
    module procedure cmlAddAtomsSP
    module procedure cmlAddAtomsSP_sh
    module procedure cmlAddAtoms_3_SP
    module procedure cmlAddAtoms_3_SP_sh
    module procedure cmlAddAtomsDP
    module procedure cmlAddAtomsDP_sh
    module procedure cmlAddAtoms_3_DP
    module procedure cmlAddAtoms_3_DP_sh
  end interface

  interface cmlAddParticles
    module procedure cmlAddParticlesSP
    module procedure cmlAddParticlesSP_sh
    module procedure cmlAddParticles_3_SP
    module procedure cmlAddParticles_3_SP_sh
    module procedure cmlAddParticlesDP
    module procedure cmlAddParticlesDP_sh
    module procedure cmlAddParticles_3_DP
    module procedure cmlAddParticles_3_DP_sh
  end interface

#ifndef DUMMYLIB
  interface cmlAddCoords
    module procedure cmlAddCoords_sp
    module procedure cmlAddCoords_dp
  end interface

  interface addDlpolyMatrix
    module procedure addDlpolyMatrix_sp
    module procedure addDlpolyMatrix_3_sp
    module procedure addDlpolyMatrix_dp
    module procedure addDlpolyMatrix_3_dp
  end interface
#endif

  public :: cmlStartMolecule
  public :: cmlEndMolecule

  public :: cmlAddAtoms
  public :: cmlAddParticles

  public :: cmlAddMolecule

contains

  subroutine cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    type(xmlf_t), intent(inout) :: xf

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



#ifndef DUMMYLIB
    
    call xml_NewElement(xf, "molecule")

    if (present(dictRef)) call xml_addAttribute(xf, "dictRef", dictRef)
    if (present(convention)) call xml_addAttribute(xf, "convention", convention)
    if (present(title)) call xml_addAttribute(xf, "title", title)
    if (present(id)) call xml_addAttribute(xf, "id", id)
    if (present(ref)) call xml_addAttribute(xf, "ref", ref)
    if (present(formula)) call xml_addAttribute(xf, "formula", formula)
    if (present(chirality)) call xml_addAttribute(xf, "chirality", chirality)
    if (present(role)) call xml_addAttribute(xf, "role", role)



#endif
  end subroutine cmlStartMolecule

  subroutine cmlEndMolecule(xf)
    type(xmlf_t), intent(inout) :: xf

#ifndef DUMMYLIB
    call xml_EndElement(xf, "molecule")
#endif
  end subroutine cmlEndMolecule

  subroutine cmlAddMoleculesp(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in)              :: coords(:, :)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    real(kind=sp), intent(in), optional :: occupancies(:)
    character(len=*), intent(in), optional :: atomIds(:)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional          :: nobondcheck

#ifndef DUMMYLIB

    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)

#endif

  end subroutine cmlAddMoleculesp


  subroutine cmlAddMoleculesp_sh(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=sp), intent(in)              :: coords(3, natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms) 
    real(kind=sp), intent(in), optional :: occupancies(natoms) 
    character(len=*), intent(in), optional :: atomIds(natoms) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional :: nobondcheck

#ifndef DUMMYLIB
    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
        call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)
#endif

  end subroutine cmlAddMoleculesp_sh


  subroutine cmlAddMolecule_3_sp(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in)              :: x(:)
    real(kind=sp), intent(in)              :: y(:)
    real(kind=sp), intent(in)              :: z(:)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    character(len=*), intent(in), optional :: atomIds(:) 
    real(kind=sp), intent(in), optional :: occupancies(:) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional :: nobondcheck

#ifndef DUMMYLIB
    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
        call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)
#endif

  end subroutine cmlAddMolecule_3_sp


  subroutine cmlAddMolecule_3_sp_sh(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=sp), intent(in)              :: x(natoms)
    real(kind=sp), intent(in)              :: y(natoms)
    real(kind=sp), intent(in)              :: z(natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms)
    character(len=*), intent(in), optional :: atomIds(natoms)
    real(kind=sp), intent(in), optional :: occupancies(natoms)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional :: nobondcheck

#ifndef DUMMYLIB
    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)
#endif

  end subroutine cmlAddMolecule_3_sp_sh

  subroutine cmlAddAtomssp(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in)              :: coords(:, :)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    real(kind=sp), intent(in), optional :: occupancies(:)
    character(len=*), intent(in), optional :: atomIds(:)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(coords,2)

    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddAtomssp

  subroutine cmlAddAtomssp_sh(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=sp), intent(in)              :: coords(3, natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms) 
    real(kind=sp), intent(in), optional :: occupancies(natoms) 
    character(len=*), intent(in), optional :: atomIds(natoms) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords(:,:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddAtomssp_sh

  subroutine cmlAddAtoms_3_sp(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in)              :: x(:)
    real(kind=sp), intent(in)              :: y(:)
    real(kind=sp), intent(in)              :: z(:)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    character(len=*), intent(in), optional :: atomIds(:) 
    real(kind=sp), intent(in), optional :: occupancies(:) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x, y, z, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(x)

    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")

#endif

  end subroutine cmlAddAtoms_3_sp

  subroutine cmlAddAtoms_3_sp_sh(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=sp), intent(in)              :: x(natoms)
    real(kind=sp), intent(in)              :: y(natoms)
    real(kind=sp), intent(in)              :: z(natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms)
    character(len=*), intent(in), optional :: atomIds(natoms)
    real(kind=sp), intent(in), optional :: occupancies(natoms)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x(:natoms), y(:natoms), z(:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddAtoms_3_sp_sh

  subroutine cmlAddParticlessp(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in)              :: coords(:, :)
    character(len=*), intent(in), optional :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    real(kind=sp), intent(in), optional :: occupancies(:)
    character(len=*), intent(in), optional :: atomIds(:)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(coords,2)

    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddParticlessp

  subroutine cmlAddParticlessp_sh(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=sp), intent(in)              :: coords(3, natoms)
    character(len=*), intent(in), optional :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms) 
    real(kind=sp), intent(in), optional :: occupancies(natoms) 
    character(len=*), intent(in), optional :: atomIds(natoms) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords(:,:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddParticlessp_sh

  subroutine cmlAddParticles_3_sp(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in)              :: x(:)
    real(kind=sp), intent(in)              :: y(:)
    real(kind=sp), intent(in)              :: z(:)
    character(len=*), intent(in), optional :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    character(len=*), intent(in), optional :: atomIds(:) 
    real(kind=sp), intent(in), optional :: occupancies(:) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x, y, z, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(x)

    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")

#endif

  end subroutine cmlAddParticles_3_sp

  subroutine cmlAddParticles_3_sp_sh(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=sp), intent(in)              :: x(natoms)
    real(kind=sp), intent(in)              :: y(natoms)
    real(kind=sp), intent(in)              :: z(natoms)
    character(len=*), intent(in), optional :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms)
    character(len=*), intent(in), optional :: atomIds(natoms)
    real(kind=sp), intent(in), optional :: occupancies(natoms)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x(:natoms), y(:natoms), z(:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddParticles_3_sp_sh

#ifndef DUMMYLIB

  subroutine cmlAddCoords_sp(xf, coords, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=sp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional  :: style
    character(len=*), intent(in), optional  :: fmt

    if (present(style)) then
      select case(style)
      case ("x3") 
        call addcoords_x3_sp(xf, coords, fmt)
      case ("cartesian") 
        call addcoords_x3_sp(xf, coords, fmt)
      case ("xFrac")
        call addcoords_xfrac_sp(xf, coords, fmt)
      case ("fractional")
        call addcoords_xfrac_sp(xf, coords, fmt)
      case ("xyz3")
        call addcoords_xyz3_sp(xf, coords, fmt)
      case ("xyzFrac")
        call addcoords_xyzfrac_sp(xf, coords, fmt)
      case default
        call FoX_error("Invalid style specification for atomic coordinates")
      end select
    else
      call addcoords_x3_sp(xf, coords, fmt)
    endif

  end subroutine cmlAddCoords_sp

  subroutine addcoords_xyz3_sp(xf, coords, fmt)
    type(xmlf_t), intent(inout)              :: xf
    real(kind=sp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional   :: fmt

    select case (size(coords))
    case (2)
      call xml_AddAttribute(xf, "xy2", coords,fmt)
    case(3)
      call xml_AddAttribute(xf, "xyz3", coords,fmt)
    end select

  end subroutine addcoords_xyz3_sp


  subroutine addcoords_xyzfrac_sp(xf, coords, fmt)
    type(xmlf_t), intent(inout)             :: xf
    real(kind=sp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional  :: fmt

    select case (size(coords))
    case (2)
      call xml_AddAttribute(xf, "xyFract", coords, fmt)
    case(3)
      call xml_AddAttribute(xf, "xyzFract", coords, fmt)
    end select

  end subroutine addcoords_xyzfrac_sp


  subroutine addcoords_x3_sp(xf, coords, fmt)
    type(xmlf_t), intent(inout)            :: xf
    real(kind=sp), intent(in), dimension(:):: coords
    character(len=*), intent(in), optional :: fmt

    select case(size(coords))
    case(2)
      call xml_AddAttribute(xf, "x2", coords(1), fmt)
      call xml_AddAttribute(xf, "y2", coords(2), fmt)
    case(3)
      call xml_AddAttribute(xf, "x3", coords(1), fmt)
      call xml_AddAttribute(xf, "y3", coords(2), fmt)
      call xml_AddAttribute(xf, "z3", coords(3), fmt)
    end select

  end subroutine addcoords_x3_sp


  subroutine addcoords_xfrac_sp(xf, coords, fmt)
    type(xmlf_t), intent(inout)             :: xf
    real(kind=sp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional  :: fmt

    call xml_AddAttribute(xf, "xFract", coords(1), fmt)
    call xml_AddAttribute(xf, "yFract", coords(2), fmt)
    call xml_AddAttribute(xf, "zFract", coords(3), fmt)

  end subroutine addcoords_xfrac_sp

  subroutine addDlpolyMatrix_sp(xf, coords, elems)
    type(xmlf_t), intent(inout)                :: xf
    real(kind=sp), intent(in), dimension(:, :) :: coords
    character(len=2), intent(in), dimension(:) :: elems

    integer :: natoms, i

    natoms = size(elems)
    
    call xml_NewElement(xf, "matrix")
    call xml_AddAttribute(xf, "nrows", size(elems))
    call xml_AddAttribute(xf, "ncols", 11)
    call xml_AddAttribute(xf, "dataType", "xsd:string")
    
    call xml_AddNewline(xf)
    do i = 1, natoms
      call xml_AddCharacters(xf, elems(i)//"  "//str(i))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, str(coords(1,i))//" "//str(coords(2,i))//" "//str(coords(3,i)))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
    enddo
    call xml_EndElement(xf, "matrix")
  end subroutine addDlpolyMatrix_sp

  subroutine addDlpolyMatrix_3_sp(xf, x, y, z, elems)
    type(xmlf_t), intent(inout)                :: xf
    real(kind=sp), intent(in), dimension(:)    :: x, y, z
    character(len=2), intent(in), dimension(:) :: elems

    integer :: natoms, i

    natoms = size(elems)
    
    call xml_NewElement(xf, "matrix")
    call xml_AddAttribute(xf, "nrows", size(elems))
    call xml_AddAttribute(xf, "ncols", 11)
    call xml_AddAttribute(xf, "dataType", "xsd:string")
    
    call xml_AddNewline(xf)
    do i = 1, natoms
      call xml_AddCharacters(xf, elems(i)//"  "//str(i))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, str(x(i))//" "//str(y(i))//" "//str(z(i)))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
    enddo
    call xml_EndElement(xf, "matrix")
  end subroutine addDlpolyMatrix_3_sp

#endif



  subroutine cmlAddMoleculedp(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in)              :: coords(:, :)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    real(kind=dp), intent(in), optional :: occupancies(:)
    character(len=*), intent(in), optional :: atomIds(:)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional          :: nobondcheck

#ifndef DUMMYLIB

    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)

#endif

  end subroutine cmlAddMoleculedp


  subroutine cmlAddMoleculedp_sh(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=dp), intent(in)              :: coords(3, natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms) 
    real(kind=dp), intent(in), optional :: occupancies(natoms) 
    character(len=*), intent(in), optional :: atomIds(natoms) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional :: nobondcheck

#ifndef DUMMYLIB
    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
        call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)
#endif

  end subroutine cmlAddMoleculedp_sh


  subroutine cmlAddMolecule_3_dp(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in)              :: x(:)
    real(kind=dp), intent(in)              :: y(:)
    real(kind=dp), intent(in)              :: z(:)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    character(len=*), intent(in), optional :: atomIds(:) 
    real(kind=dp), intent(in), optional :: occupancies(:) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional :: nobondcheck

#ifndef DUMMYLIB
    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
        call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)
#endif

  end subroutine cmlAddMolecule_3_dp


  subroutine cmlAddMolecule_3_dp_sh(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt &
,dictRef,convention,title,id,ref,formula,chirality,role , &
bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds, nobondcheck)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=dp), intent(in)              :: x(natoms)
    real(kind=dp), intent(in)              :: y(natoms)
    real(kind=dp), intent(in)              :: z(natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms)
    character(len=*), intent(in), optional :: atomIds(natoms)
    real(kind=dp), intent(in), optional :: occupancies(natoms)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

    character(len=*), intent(in), optional :: dictRef
    character(len=*), intent(in), optional :: convention
    character(len=*), intent(in), optional :: title
    character(len=*), intent(in), optional :: id
    character(len=*), intent(in), optional :: ref
    character(len=*), intent(in), optional :: formula
    character(len=*), intent(in), optional :: chirality
    character(len=*), intent(in), optional :: role



    character(len=*), intent(in), optional :: bondAtom1Refs(:)
    character(len=*), intent(in), optional :: bondAtom2Refs(:)
    character(len=*), intent(in), optional :: bondOrders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    logical, intent(in), optional :: nobondcheck

#ifndef DUMMYLIB
    call cmlStartMolecule(xf &
,dictRef,convention,title,id,ref,formula,chirality,role)
    call cmlAddAtoms(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    if (present(bondAtom1Refs)) then
      if (present(bondAtom2Refs).and.present(bondOrders)) then
        if (present(atomIds)) then
          call checkBondIdRefs(atomIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
          call addBondArray(xf, bondAtom1Refs, bondAtom2Refs, bondOrders, bondIds)
        else
          call FoX_error("AtomIds must be provided to add bonds")
        endif
      else
        call FoX_error("Two AtomRefs arrays and a bondOrder array must be provided to add bonds")
      endif
    endif
    call cmlEndMolecule(xf)
#endif

  end subroutine cmlAddMolecule_3_dp_sh

  subroutine cmlAddAtomsdp(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in)              :: coords(:, :)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    real(kind=dp), intent(in), optional :: occupancies(:)
    character(len=*), intent(in), optional :: atomIds(:)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(coords,2)

    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddAtomsdp

  subroutine cmlAddAtomsdp_sh(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=dp), intent(in)              :: coords(3, natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms) 
    real(kind=dp), intent(in), optional :: occupancies(natoms) 
    character(len=*), intent(in), optional :: atomIds(natoms) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords(:,:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddAtomsdp_sh

  subroutine cmlAddAtoms_3_dp(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in)              :: x(:)
    real(kind=dp), intent(in)              :: y(:)
    real(kind=dp), intent(in)              :: z(:)
    character(len=*), intent(in)           :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    character(len=*), intent(in), optional :: atomIds(:) 
    real(kind=dp), intent(in), optional :: occupancies(:) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x, y, z, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(x)

    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")

#endif

  end subroutine cmlAddAtoms_3_dp

  subroutine cmlAddAtoms_3_dp_sh(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=dp), intent(in)              :: x(natoms)
    real(kind=dp), intent(in)              :: y(natoms)
    real(kind=dp), intent(in)              :: z(natoms)
    character(len=*), intent(in)           :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms)
    character(len=*), intent(in), optional :: atomIds(natoms)
    real(kind=dp), intent(in), optional :: occupancies(natoms)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x(:natoms), y(:natoms), z(:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "atom")
      call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "atom")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddAtoms_3_dp_sh

  subroutine cmlAddParticlesdp(xf, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in)              :: coords(:, :)
    character(len=*), intent(in), optional :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    real(kind=dp), intent(in), optional :: occupancies(:)
    character(len=*), intent(in), optional :: atomIds(:)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(coords,2)

    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddParticlesdp

  subroutine cmlAddParticlesdp_sh(xf, natoms, elements, atomRefs, coords, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=dp), intent(in)              :: coords(3, natoms)
    character(len=*), intent(in), optional :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms) 
    real(kind=dp), intent(in), optional :: occupancies(natoms) 
    character(len=*), intent(in), optional :: atomIds(natoms) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, coords(:,:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=coords(:,i), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddParticlesdp_sh

  subroutine cmlAddParticles_3_dp(xf, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in)              :: x(:)
    real(kind=dp), intent(in)              :: y(:)
    real(kind=dp), intent(in)              :: z(:)
    character(len=*), intent(in), optional :: elements(:)
    character(len=*), intent(in), optional :: atomRefs(:) 
    character(len=*), intent(in), optional :: atomIds(:) 
    real(kind=dp), intent(in), optional :: occupancies(:) 
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i, natoms

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x, y, z, elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")

    natoms = size(x)

    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")

#endif

  end subroutine cmlAddParticles_3_dp

  subroutine cmlAddParticles_3_dp_sh(xf, natoms, elements, x, y, z, atomRefs, occupancies, atomIds, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    integer, intent(in) :: natoms
    real(kind=dp), intent(in)              :: x(natoms)
    real(kind=dp), intent(in)              :: y(natoms)
    real(kind=dp), intent(in)              :: z(natoms)
    character(len=*), intent(in), optional :: elements(natoms)
    character(len=*), intent(in), optional :: atomRefs(natoms)
    character(len=*), intent(in), optional :: atomIds(natoms)
    real(kind=dp), intent(in), optional :: occupancies(natoms)
    character(len=*), intent(in), optional :: fmt
    character(len=*), intent(in), optional :: style

#ifndef DUMMYLIB
    integer          :: i

    if (present(style)) then
      if (style=="DL_POLY") then
        if (present(atomRefs).or.present(occupancies).or.present(atomIds).or.present(fmt)) &
          call FoX_error("With DL_POLY style, no optional arguments permitted.")
        call addDlpolyMatrix(xf, x(:natoms), y(:natoms), z(:natoms), elements)
        return
      endif
    endif


    call xml_NewElement(xf, "atomArray")


    do i = 1, natoms
      call xml_NewElement(xf, "particle")
      if (present(elements)) call xml_AddAttribute(xf, "elementType", trim(elements(i)))
      call cmlAddCoords(xf, coords=(/x(i),y(i),z(i)/), style=style, fmt=fmt)
      if (present(occupancies)) call xml_AddAttribute(xf, "occupancy", occupancies(i))
      if (present(atomRefs)) call xml_AddAttribute(xf, "ref", trim(atomRefs(i)))
      if (present(atomIds)) call xml_AddAttribute(xf, "id", trim(atomIds(i)))
      call xml_EndElement(xf, "particle")
     enddo


    call xml_EndElement(xf, "atomArray")
#endif

  end subroutine cmlAddParticles_3_dp_sh

#ifndef DUMMYLIB

  subroutine cmlAddCoords_dp(xf, coords, style, fmt)
    type(xmlf_t), intent(inout) :: xf
    real(kind=dp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional  :: style
    character(len=*), intent(in), optional  :: fmt

    if (present(style)) then
      select case(style)
      case ("x3") 
        call addcoords_x3_dp(xf, coords, fmt)
      case ("cartesian") 
        call addcoords_x3_dp(xf, coords, fmt)
      case ("xFrac")
        call addcoords_xfrac_dp(xf, coords, fmt)
      case ("fractional")
        call addcoords_xfrac_dp(xf, coords, fmt)
      case ("xyz3")
        call addcoords_xyz3_dp(xf, coords, fmt)
      case ("xyzFrac")
        call addcoords_xyzfrac_dp(xf, coords, fmt)
      case default
        call FoX_error("Invalid style specification for atomic coordinates")
      end select
    else
      call addcoords_x3_dp(xf, coords, fmt)
    endif

  end subroutine cmlAddCoords_dp

  subroutine addcoords_xyz3_dp(xf, coords, fmt)
    type(xmlf_t), intent(inout)              :: xf
    real(kind=dp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional   :: fmt

    select case (size(coords))
    case (2)
      call xml_AddAttribute(xf, "xy2", coords,fmt)
    case(3)
      call xml_AddAttribute(xf, "xyz3", coords,fmt)
    end select

  end subroutine addcoords_xyz3_dp


  subroutine addcoords_xyzfrac_dp(xf, coords, fmt)
    type(xmlf_t), intent(inout)             :: xf
    real(kind=dp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional  :: fmt

    select case (size(coords))
    case (2)
      call xml_AddAttribute(xf, "xyFract", coords, fmt)
    case(3)
      call xml_AddAttribute(xf, "xyzFract", coords, fmt)
    end select

  end subroutine addcoords_xyzfrac_dp


  subroutine addcoords_x3_dp(xf, coords, fmt)
    type(xmlf_t), intent(inout)            :: xf
    real(kind=dp), intent(in), dimension(:):: coords
    character(len=*), intent(in), optional :: fmt

    select case(size(coords))
    case(2)
      call xml_AddAttribute(xf, "x2", coords(1), fmt)
      call xml_AddAttribute(xf, "y2", coords(2), fmt)
    case(3)
      call xml_AddAttribute(xf, "x3", coords(1), fmt)
      call xml_AddAttribute(xf, "y3", coords(2), fmt)
      call xml_AddAttribute(xf, "z3", coords(3), fmt)
    end select

  end subroutine addcoords_x3_dp


  subroutine addcoords_xfrac_dp(xf, coords, fmt)
    type(xmlf_t), intent(inout)             :: xf
    real(kind=dp), intent(in), dimension(:) :: coords
    character(len=*), intent(in), optional  :: fmt

    call xml_AddAttribute(xf, "xFract", coords(1), fmt)
    call xml_AddAttribute(xf, "yFract", coords(2), fmt)
    call xml_AddAttribute(xf, "zFract", coords(3), fmt)

  end subroutine addcoords_xfrac_dp

  subroutine addDlpolyMatrix_dp(xf, coords, elems)
    type(xmlf_t), intent(inout)                :: xf
    real(kind=dp), intent(in), dimension(:, :) :: coords
    character(len=2), intent(in), dimension(:) :: elems

    integer :: natoms, i

    natoms = size(elems)
    
    call xml_NewElement(xf, "matrix")
    call xml_AddAttribute(xf, "nrows", size(elems))
    call xml_AddAttribute(xf, "ncols", 11)
    call xml_AddAttribute(xf, "dataType", "xsd:string")
    
    call xml_AddNewline(xf)
    do i = 1, natoms
      call xml_AddCharacters(xf, elems(i)//"  "//str(i))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, str(coords(1,i))//" "//str(coords(2,i))//" "//str(coords(3,i)))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
    enddo
    call xml_EndElement(xf, "matrix")
  end subroutine addDlpolyMatrix_dp

  subroutine addDlpolyMatrix_3_dp(xf, x, y, z, elems)
    type(xmlf_t), intent(inout)                :: xf
    real(kind=dp), intent(in), dimension(:)    :: x, y, z
    character(len=2), intent(in), dimension(:) :: elems

    integer :: natoms, i

    natoms = size(elems)
    
    call xml_NewElement(xf, "matrix")
    call xml_AddAttribute(xf, "nrows", size(elems))
    call xml_AddAttribute(xf, "ncols", 11)
    call xml_AddAttribute(xf, "dataType", "xsd:string")
    
    call xml_AddNewline(xf)
    do i = 1, natoms
      call xml_AddCharacters(xf, elems(i)//"  "//str(i))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, str(x(i))//" "//str(y(i))//" "//str(z(i)))
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
      call xml_AddCharacters(xf, "0 0 0")
      call xml_AddNewline(xf)
    enddo
    call xml_EndElement(xf, "matrix")
  end subroutine addDlpolyMatrix_3_dp

#endif



#ifndef DUMMYLIB
  subroutine addBondArray(xf, atom1Refs, atom2Refs, orders, bondIds)
    type(xmlf_t), intent(inout)            :: xf
    character(len=*), intent(in)           :: atom1Refs(:)
    character(len=*), intent(in)           :: atom2Refs(:)
    character(len=*), intent(in)           :: orders(:)
    character(len=*), intent(in), optional :: bondIds(:)
    integer                                :: nbonds
    integer                                :: i

    nbonds = size(atom1Refs)
    ! Basic argument verification
    if (size(atom2Refs).ne.nbonds) &
      call FoX_error("Length of atomRef arrays must match in WCML addBondArray")
    if (size(orders).ne.nbonds) &
      call FoX_error("Length of atomRef and order arrays must match in WCML addBondArray")
    if (present(bondIds)) then
      if (size(bondIds).ne.nbonds) &
        call FoX_error("Length of atomRef and bondId arrays must match in WCML addBondArray")
    endif

    ! Add the bond array
    call xml_NewElement(xf, "bondArray")
       do i = 1, nbonds
         call xml_NewElement(xf, "bond")
             call xml_AddAttribute(xf, "atomRefs2", trim(atom1Refs(i))//" "//trim(atom2Refs(i)))
             call xml_AddAttribute(xf, "order", orders(i))
             if (present(bondIds)) & 
               call xml_AddAttribute(xf, "id", trim(bondIds(i))) 
         call xml_EndElement(xf, "bond")
       enddo 
    call xml_EndElement(xf, "bondArray")

  end subroutine addBondArray

  subroutine checkBondIdRefs(atomArrayIds, bondAtom1Refs, bondAtom2Refs, nobondcheck)
    character(len=*), intent(in)           :: atomArrayIds(:)
    character(len=*), intent(in)           :: bondAtom1Refs(:)
    character(len=*), intent(in)           :: bondAtom2Refs(:)
    logical, intent(in), optional          :: nobondcheck
    logical                                :: bondmatrix(size(atomArrayIds),size(atomArrayIds))
    integer                                :: nbonds
    integer                                :: natoms
    integer                                :: i
    integer                                :: j
    logical                                :: bond1OK
    logical                                :: bond2OK
    integer                                :: atom1num
    integer                                :: atom2num

    if (present(nobondcheck)) then
      if (nobondcheck) return ! skip all checks
    endif

    atom1num = -1 ! Supress bogus gfortran warning
    atom2num = -1
    bondmatrix = .false.
    natoms = size(atomArrayIds)
    nbonds = size(bondAtom1Refs)
    if (size(bondAtom2Refs).ne.nbonds) &
      call FoX_error("Length of atomRef arrays must match in WCML checkBondIdRefs")

    do i = 1, nbonds
      if (trim(bondAtom1Refs(i)).eq.trim(bondAtom2Refs(i))) &
        call FoX_error("The two atomRefs in a bond must be different")
      bond1OK = .false.
      bond2OK = .false.
      do j = 1, natoms
        if (trim(bondAtom1Refs(i)).eq.trim(atomArrayIds(j))) &
          bond1OK = .true.
          atom1num = j
        if (trim(bondAtom2Refs(i)).eq.trim(atomArrayIds(j))) &
          bond2OK = .true.
          atom2num = j
        if (bond1OK.and.bond2OK) exit
      enddo
      if (.not.bond1OK) call FoX_error(bondAtom1Refs(i) // " not found in checkBondIdRefs")
      if (.not.bond2OK) call FoX_error(bondAtom2Refs(i) // " not found in checkBondIdRefs")
      ! Both atoms bust have been found to get here...
      if (bondmatrix(atom1num,atom2num)) then
        ! Seen this bond before
        call FoX_error("A bond cannot be added twice.")
      else
        ! We've seen this bond (both ways) - so don't forget
        bondmatrix(atom1num,atom2num) = .true.
        bondmatrix(atom2num,atom1num) = .true.
      endif
    enddo

  end subroutine checkBondIdRefs
#endif

end module m_wcml_molecule
