anchor.diff:       & plotdef%box%origin%coord
         v3 (:) = plotdef%box%pointarray(3)%point%coord - &
        & plotdef%box%origin%coord
  +
         ip = 0
  -      Do ip1 = 0, plotdef%box%grid(1) - 1
  -         t1 = dble (ip1) / dble (plotdef%box%grid(1))
  +      Do ip3 = 0, plotdef%box%grid(3) - 1
  +         t3 = dble (ip3) / dble (plotdef%box%grid(3))
            Do ip2 = 0, plotdef%box%grid(2) - 1
               t2 = dble (ip2) / dble (plotdef%box%grid(2))
  -            Do ip3 = 0, plotdef%box%grid(3) - 1
  -               t3 = dble (ip3) / dble (plotdef%box%grid(3))
  +            Do ip1 = 0, plotdef%box%grid(1) - 1
  +               t1 = dble (ip1) / dble (plotdef%box%grid(1))
                  ip = ip + 1
                  vpl (:, ip) = t1 * v1 (:) + t2 * v2 (:) + t3 * v3 (:) + &
                 & plotdef%box%origin%coord
  
positive.code:         & fp(:, i))
        End Do
  !write xml
     Call xml_OpenFile (fname//"3d.xml", xf, replace=.True., &
       & pretty_print=.True.)
        Call xml_NewElement (xf, "plot3d")
        Call xml_NewElement (xf, "title")
  
negative.code:       Call xml_NewElement (xf, "grid")
        Write (buffer, '(3I6)') plotdef%box%grid
        Call xml_AddAttribute (xf, "gridticks", trim(adjustl(buffer)))
        Write (buffer, '(6G18.10)') plotdef%box%origin%coord
        Call xml_AddAttribute (xf, "origin", trim(adjustl(buffer)))
        !write x axis description
         call DGEMV('N',3,3,1.d0, input%structure%crystal%basevect(1,1),3,&
        plotdef%box%origin%coord,1,0.d0,tmpv(1),1)
        Write (buffer, '(3G18.10)') tmpv
        Call xml_AddAttribute (xf, "originrs", trim(adjustl(buffer)))
        Call xml_NewElement (xf, "axis")
        Call xml_AddAttribute (xf, "name", "a")
        Write (buffer, '(6G18.10)') plotdef%box%pointarray(1)%point%coord
        Call xml_AddAttribute (xf, "endpoint", trim(adjustl(buffer)))
        Write (buffer, '(6G18.10)') &
       & (plotdef%box%pointarray(1)%point%coord-plotdef%box%origin%coord) &
       & / plotdef%box%grid(1)
        Call xml_AddAttribute (xf, "delta", trim(adjustl(buffer)))
        call DGEMV('N',3,3,1.d0, input%structure%crystal%basevect(1,1),3,&
        plotdef%box%pointarray(1)%point%coord,1,0.d0,tmpv(1),1)
        Write (buffer, '(3G18.10)') tmpv
        Call xml_AddAttribute (xf, "endpointrs", trim(adjustl(buffer)))
        Call xml_endElement (xf, "axis")
        !write y axis description
        Call xml_NewElement (xf, "axis")
        Call xml_AddAttribute (xf, "name", "b")
        Write (buffer, '(6G18.10)') plotdef%box%pointarray(2)%point%coord
        Call xml_AddAttribute (xf, "endpoint", trim(adjustl(buffer)))
        Write (buffer, '(6G18.10)') &
       & (plotdef%box%pointarray(2)%point%coord-plotdef%box%origin%coord) &
       & / plotdef%box%grid(2)
        Call xml_AddAttribute (xf, "delta", trim(adjustl(buffer)))
  
        call DGEMV('N',3,3,1.d0, input%structure%crystal%basevect(1,1),3,&
        plotdef%box%pointarray(2)%point%coord,1,0.d0,tmpv(1),1)
        Write (buffer, '(3G18.10)') tmpv
        Call xml_AddAttribute (xf, "endpointrs", trim(adjustl(buffer)))
  !
        Call xml_endElement (xf, "axis")
              !write z axis description
        Call xml_NewElement (xf, "axis")
        Call xml_AddAttribute (xf, "name", "c")
        Write (buffer, '(6G18.10)') plotdef%box%pointarray(3)%point%coord
        Call xml_AddAttribute (xf, "endpoint", trim(adjustl(buffer)))
        Write (buffer, '(6G18.10)') &
       & (plotdef%box%pointarray(3)%point%coord-plotdef%box%origin%coord) &
       & / plotdef%box%grid(3)
        Call xml_AddAttribute (xf, "delta", trim(adjustl(buffer)))
         call DGEMV('N',3,3,1.d0, input%structure%crystal%basevect(1,1),3,&
        plotdef%box%pointarray(3)%point%coord,1,0.d0,tmpv(1),1)
        Write (buffer, '(3G18.10)') tmpv
        Call xml_AddAttribute (xf, "endpointrs", trim(adjustl(buffer)))
  !
        Call xml_endElement (xf, "axis")
        Call xml_endElement (xf, "grid")
  !
  !
  ! write functions to file
        ip = 0
        Do ifunction = 1, nf
           Call xml_NewElement (xf, "function")
           Write (buffer20, '(I14)') np
           Call xml_AddAttribute (xf, "n", trim(adjustl(buffer20)))
           Do ip1 = 0, plotdef%box%grid(3) - 1
              Call xml_NewElement (xf, "row")
              Call xml_AddAttribute (xf, "const", "c")
              Write (buffer20, '(I14)') ip1
              Call xml_AddAttribute (xf, "index", &
             & trim(adjustl(buffer20)))
               Do ip2 = 0, plotdef%box%grid(2) - 1
                 Call xml_NewElement (xf, "row")
                 Call xml_AddAttribute (xf, "const", "b")
                 Write (buffer20, '(I14)') ip2
                 Call xml_AddAttribute (xf, "index", &
                & trim(adjustl(buffer20)))
                 Do ip3 = 0, plotdef%box%grid(1) - 1
                    ip = ip + 1
                    Write (buffer20, '(6G18.10)') fp (ip, ifunction)
                    Call xml_AddCharacters (xf, buffer20)
                 End Do
                 Call xml_endElement (xf, "row")
              End Do
              Call xml_endElement (xf, "row")
           End Do
  !
           Call xml_NewElement (xf, "function")
        End Do
  !
  !
        Deallocate (vpl, fp)
        Call xml_Close (xf)
  