text:
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
! H0 X
! H0 X   libAtoms+QUIP: atomistic simulation library
! H0 X
! H0 X   Portions of this code were written by
! H0 X     Albert Bartok-Partay, Silvia Cereda, Gabor Csanyi, James Kermode,
! H0 X     Ivan Solt, Wojciech Szlachta, Csilla Varnai, Steven Winfield.
! H0 X
! H0 X   Copyright 2006-2010.
! H0 X
! H0 X   These portions of the source code are released under the GNU General
! H0 X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html
! H0 X
! H0 X   If you would like to license the source code under different terms,
! H0 X   please contact Gabor Csanyi, gabor@csanyi.net
! H0 X
! H0 X   Portions of this code were written by Noam Bernstein as part of
! H0 X   his employment for the U.S. Government, and are not subject
! H0 X   to copyright in the USA.
! H0 X
! H0 X
! H0 X   When using this software, please cite the following reference:
! H0 X
! H0 X   http://www.libatoms.org
! H0 X
! H0 X  Additional contributions by
! H0 X    Alessio Comisso, Chiara Gattinoni, and Gianpietro Moras
! H0 X
! H0 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

!X
!X  IPModel_SW_VP  module
!X
!%  Danny potential J ChemPhys 127, 204704(2007) for SiO_2.
!%
!%  Only the case of fourfold coordinated silcon and twofold coordinated
!%  is considered. (-->Fixed charges).
!%
!%  The IPModel_SW_VP object contains all the parameters read
!%  from an 'SW_VP_params' XML stanza.
!X
!X  Contribution from Anke Butenuth
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
#include "error.inc"

module IPModel_SW_VP_module

use error_module
use system_module, only : dp, inoutput, print, PRINT_NERD, PRINT_ANALYSIS, current_verbosity, verbosity_push_decrement, verbosity_pop, operator(//)
use dictionary_module
use paramreader_module
use linearalgebra_module
use atoms_types_module
use atoms_module

use mpi_context_module
use QUIP_Common_module
use Functions_module
implicit none

private

include 'IPModel_interface.h'

public :: IPModel_SW_VP
type IPModel_SW_VP
  integer :: n_types = 0         !% Number of atomic types
  integer, allocatable :: atomic_num(:), type_of_atomic_num(:)  !% Atomic number dimensioned as \texttt{n_types}

  real(dp) :: cutoff = 0.0_dp

  real(dp), allocatable :: a(:,:), AA(:,:), BB(:,:), p(:,:), q(:,:), sigma(:,:), eps2(:,:), C_0(:,:), C_1(:,:), b(:,:), D_SiO(:,:), C_OOprime(:,:), D_OOprime(:,:) !% IP parameters
  real(dp), allocatable :: lambda(:,:,:), d1(:,:,:), d2(:,:,:), gamma1(:,:,:), gamma2(:,:,:), eps3(:,:,:), costheta0(:,:,:) !% IP parameters

  character(len=STRING_LENGTH) :: label

end type IPModel_SW_VP

logical, private :: parse_in_ip, parse_matched_label
type(IPModel_SW_VP), private, pointer :: parse_ip

interface Initialise
  module procedure IPModel_SW_VP_Initialise_str
end interface Initialise

interface Finalise
  module procedure IPModel_SW_VP_Finalise
end interface Finalise

interface Print
  module procedure IPModel_SW_VP_Print
end interface Print

interface Calc
  module procedure IPModel_SW_VP_Calc
end interface Calc

contains

subroutine IPModel_SW_VP_Initialise_str(this, args_str, param_str)
  type(IPModel_SW_VP), intent(inout) :: this
  character(len=*), intent(in) :: args_str, param_str

  type(Dictionary) :: params

  call Finalise(this)

  call initialise(params)
  this%label=''
  call param_register(params, 'label', '', this%label, help_string="No help yet.  This source file was $LastChangedBy: nb326 $")
  if (.not. param_read_line(params, args_str, ignore_unknown=.true.,task='IPModel_SW_VP_Initialise_str args_str')) then
    call system_abort("IPModel_SW_VP_Initialise_str failed to parse label from args_str="//trim(args_str))
  endif
  call finalise(params)

  call IPModel_SW_VP_read_params_xml(this, param_str)

end subroutine IPModel_SW_VP_Initialise_str

subroutine IPModel_SW_VP_Finalise(this)
  type(IPModel_SW_VP), intent(inout) :: this

  if (allocated(this%atomic_num)) deallocate(this%atomic_num)
  if (allocated(this%type_of_atomic_num)) deallocate(this%type_of_atomic_num)

  if (allocated(this%a)) deallocate(this%a)
  if (allocated(this%AA)) deallocate(this%AA)
  if (allocated(this%BB)) deallocate(this%BB)
  if (allocated(this%p)) deallocate(this%p)
  if (allocated(this%q)) deallocate(this%q)
  if (allocated(this%sigma)) deallocate(this%sigma)
  if (allocated(this%eps2)) deallocate(this%eps2)
  if (allocated(this%lambda)) deallocate(this%lambda)
  if (allocated(this%gamma1)) deallocate(this%gamma1)
  if (allocated(this%gamma2)) deallocate(this%gamma2)
  if (allocated(this%eps3)) deallocate(this%eps3)
  if (allocated(this%costheta0)) deallocate(this%costheta0)
  if (allocated(this%d1)) deallocate(this%d1)
  if (allocated(this%d2)) deallocate(this%d2)
  if (allocated(this%C_0)) deallocate(this%C_0)
  if (allocated(this%C_1)) deallocate(this%C_1)
  if (allocated(this%b)) deallocate(this%b)
  if (allocated(this%D_SiO)) deallocate(this%D_SiO)
  if (allocated(this%C_OOprime)) deallocate(this%C_OOprime)
  if (allocated(this%D_OOprime)) deallocate(this%D_OOprime)

  this%n_types = 0
  this%label = ''
end subroutine IPModel_SW_VP_Finalise

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!% The potential calculator. It computes energy, forces and virial.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

subroutine IPModel_SW_VP_Calc(this, at, e, local_e, f, virial, local_virial, args_str, mpi, error)
  type(IPModel_SW_VP), intent(inout) :: this
  type(Atoms), intent(in) :: at
  real(dp), intent(out), optional :: e, local_e(:) !% \texttt{e} = System total energy, \texttt{local_e} = energy of each atom, vector dimensioned as \texttt{at%N}.
  real(dp), intent(out), optional :: f(:,:), local_virial(:,:)   !% Forces, dimensioned as \texttt{f(3,at%N)}, local virials, dimensioned as \texttt{local_virial(9,at%N)}
  real(dp), intent(out), optional :: virial(3,3)   !% Virial
  character(len=*), optional      :: args_str
  type(MPI_Context), intent(in), optional :: mpi
  integer, intent(out), optional :: error

  real(dp), pointer :: w_e(:)
  integer i, ji, j, ki, k
  real(dp) :: drij(3), drij_mag, drik(3), drik_mag, drij_dot_drik
  real(dp) :: w_f

  integer ti, tj, tk

  real(dp) :: drij_dri(3), drij_drj(3), drik_dri(3), drik_drk(3)
  real(dp) :: dcos_ijk_dri(3), dcos_ijk_drj(3), dcos_ijk_drk(3)
  real(dp) :: virial_i(3,3)

  real(dp) :: de, de_dr, de_drij, de_drik, de_dcos_ijk
  real(dp) :: cur_cutoff

  integer :: n_neigh_i

  type(Dictionary)                :: params
  logical :: has_atom_mask_name, atom_mask_exclude_all
  character(STRING_LENGTH) :: atom_mask_name
  logical, dimension(:), pointer :: atom_mask_pointer

#ifdef _OPENMP
  real(dp) :: private_virial(3,3), private_e
  real(dp), allocatable :: private_f(:,:), private_local_e(:), private_local_virial(:,:)
#endif

  INIT_ERROR(error)

  call print("IPModel_SW_VP_Calc starting ", PRINT_ANALYSIS)
  if (present(e)) e = 0.0_dp

  if (present(local_e)) then
     call check_size('Local_E',local_e,(/at%N/),'IPModel_SW_VP_Calc', error)
     local_e = 0.0_dp
  endif

  if (present(f)) then
     call check_size('Force',f,(/3,at%N/),'IPModel_SW_VP_Calc', error)
     f = 0.0_dp
  end if

  if (present(virial)) virial = 0.0_dp

  if (present(local_virial)) then
     call check_size('Local_virial',local_virial,(/9,at%N/),'IPModel_GAP_Calc', error)
     local_virial = 0.0_dp
  endif

  atom_mask_pointer => null()
  if (present(args_str)) then
     call initialise(params)
     call param_register(params, 'atom_mask_name', 'NONE', atom_mask_name, has_value_target=has_atom_mask_name, help_string="Name of logical property used to mask atoms")
     call param_register(params, 'atom_mask_exclude_all', 'F', atom_mask_exclude_all, help_string="If true, exclude contributions made by atoms with atom_mask==.false. on their neighbour with atom_mask==.true., i.e. behave as if excluded atoms were not there at all")

     if(.not. param_read_line(params, args_str, ignore_unknown=.true.,task='IPModel_SW_VP_Calc args_str')) then
        RAISE_ERROR("IPModel_SW_VP_Calc failed to parse args_str='"//trim(args_str)//"'",error)
     endif
     call finalise(params)
     if(has_atom_mask_name) then
        if (.not. assign_pointer(at, trim(atom_mask_name) , atom_mask_pointer)) &
        call system_abort("IPModel_SW_VP_Calc did not find "//trim(atom_mask_name)//" property in the atoms object.")
     else
        atom_mask_pointer => null()
     endif
  endif

  if (.not.assign_pointer(at,"weight", w_e)) nullify(w_e)

#ifdef _OPENMP
      !$omp parallel default(none) private(i, ji, j, ki, k, drij, drij_mag, drik, drik_mag, drij_dot_drik, w_f, ti, tj, tk, drij_dri, drij_drj, drik_dri, drik_drk, dcos_ijk_dri, dcos_ijk_drj, dcos_ijk_drk, de, de_dr, de_drij, de_drik, de_dcos_ijk, cur_cutoff, private_virial, private_e, private_f, private_local_e, private_local_virial, n_neigh_i, virial_i) shared(this,at,e,local_e,f,virial,local_virial,mpi,atom_mask_pointer,atom_mask_exclude_all,w_e)

  if (present(e)) private_e = 0.0_dp
  if (present(local_e)) then
    allocate(private_local_e(at%N))
    private_local_e = 0.0_dp
  endif
  if (present(f)) then
    allocate(private_f(3,at%N))
    private_f = 0.0_dp
  endif
  if (present(virial)) private_virial = 0.0_dp
  if (present(local_virial)) then
    allocate(private_local_virial(9,at%N))
    private_local_virial = 0.0_dp
  endif

!$omp do
#endif
  do i=1, at%N
   if (present(mpi)) then
       if (mpi%active) then
         if (mod(i-1, mpi%n_procs) /= mpi%my_proc) cycle
       endif
    endif

    if(associated(atom_mask_pointer)) then
       if(.not. atom_mask_pointer(i)) cycle
    endif

    ti = get_type(this%type_of_atomic_num, at%Z(i))

    if (current_verbosity() >= PRINT_ANALYSIS) call print ("IPModel_SW_VP_Calc i " // i // " " // n_neighbours(at,i), PRINT_ANALYSIS)
    cur_cutoff = maxval(this%a(ti,:)*this%sigma(ti,:))
    n_neigh_i = n_neighbours(at, i)

    do ji=1, n_neigh_i
      j = neighbour(at, i, ji, drij_mag, cosines = drij, max_dist = cur_cutoff)

      if(associated(atom_mask_pointer) .and. atom_mask_exclude_all) then
         if(.not. atom_mask_pointer(j)) cycle
      endif

      if (j <= 0) cycle
      if (drij_mag .feq. 0.0_dp) cycle

      tj = get_type(this%type_of_atomic_num, at%Z(j))

      if (drij_mag/this%sigma(ti,tj)  > this%a(ti,tj)) cycle

      if (current_verbosity() >= PRINT_ANALYSIS) call print ("IPModel_SW_VP_Calc i j " // i // " " // j, PRINT_ANALYSIS)

      if (associated(w_e)) then
        w_f = 0.5_dp*(w_e(i)+w_e(j))  !
      else
        w_f = 1.0_dp
      endif

      if (present(e) .or. present(local_e)) then
        ! factor of 0.5 because SW_VP definition goes over each pair only once

        ! select functional form of pair potentials depending on species : oxygen -oxygen
        if ((ti .eq. 2) .and. (tj  .eq. 2)) then  !
          de = 0.5_dp*this%eps2(ti,tj)*f2OO(this, drij_mag, ti, tj)

        elseif ((ti .eq. 2) .and. (tj .eq. 3)) then  ! Oxygen Silicon
          de = 0.5_dp*this%eps2(ti,tj)*f2SiO(this, drij_mag, ti, tj)

        elseif ((ti .eq. 3) .and. (tj .eq.2 )) then ! Calls the Silicon-Oxygen function again
          de = 0.5_dp*this%eps2(ti,tj)*f2SiO(this, drij_mag, ti, tj)

        elseif ((ti .eq. 3) .and. (tj .eq. 3)) then ! Silicon -Silicon
          de = 0.5_dp*this%eps2(ti,tj)*f2SiSi(this, drij_mag, ti, tj)

        else
          de= 0.0_dp
        endif

        if (present(local_e)) then
#ifdef _OPENMP
          private_local_e(i) = private_local_e(i) + de
#else
          local_e(i) = local_e(i) + de
#endif
        endif
        if (present(e)) then
#ifdef _OPENMP
          private_e = private_e + de*w_f
#else
          e = e + de*w_f
#endif
        endif
      endif

      if (present(f) .or. present(virial) .or. present(local_virial)) then

        !select pair potentials depending on species : oxygen-oxygen
        if ((ti .eq. 2) .and. (tj .eq. 2)) then
        de_dr = 0.5_dp*this%eps2(ti,tj)*df2OO_dr(this, drij_mag, ti, tj)

        elseif ((ti .eq. 3) .and. (tj .eq. 2)) then  ! oxygen-silicon
        de_dr = 0.5_dp*this%eps2(ti,tj)* df2SiO_dr(this, drij_mag, ti, tj)

        elseif ((ti.eq. 2) .and. (tj .eq. 3 )) then ! Calls the silicon-oxygen function again
        de_dr = 0.5_dp*this%eps2(ti,tj)* df2SiO_dr(this, drij_mag, ti, tj)

        elseif ((ti .eq. 3 ) .and. (tj .eq. 3)) then ! silicon-silicon
        de_dr = 0.5_dp*this%eps2(ti,tj)*df2SiSi_dr(this, drij_mag, ti, tj)
        else
        de_dr = 0.0_dp
        endif

        if (present(f)) then
#ifdef _OPENMP
          private_f(:,i) = private_f(:,i) + de_dr*w_f*drij
          private_f(:,j) = private_f(:,j) - de_dr*w_f*drij
#else
          f(:,i) = f(:,i) + de_dr*w_f*drij
          f(:,j) = f(:,j) - de_dr*w_f*drij
#endif
        endif

        if(present(virial) .or. present(local_virial)) virial_i = de_dr*w_f*(drij .outer. drij)*drij_mag

        if (present(virial)) then
#ifdef _OPENMP
          private_virial = private_virial - virial_i
#else
          virial = virial - virial_i
#endif
        endif
        if (present(local_virial)) then
#ifdef _OPENMP
           private_local_virial(:,i) = private_local_virial(:,i) - reshape(virial_i, (/9/))
#else
           local_virial(:,i) = local_virial(:,i) - reshape(virial_i, (/9/))
#endif
        endif

      endif

      if (associated(w_e)) then
        w_f = w_e(i)
      else
        w_f = 1.0_dp
      endif
      enddo !  closes j loop
      enddo !  closes i loop

      ! New loop for three body terms, because (see parameters/ip.params.SW_VP.xml) cutoff for three body terms is in SiSiO case larger then the  SiSi cutoff

      !$omp do
      do i=1, at%N

        if (present(mpi)) then
           if (mpi%active) then
             if (mod(i-1, mpi%n_procs) /= mpi%my_proc) cycle
           endif
        endif

        if(associated(atom_mask_pointer)) then
           if(.not. atom_mask_pointer(i)) cycle
        endif

    ti = get_type(this%type_of_atomic_num, at%Z(i)) !

   ! if (current_verbosity() >= PRINT_ANALYSIS) call print ("IPModel_SW_VP_Calc i " // i // " " // n_neighbours(at,i), PRINT_ANALYSIS)
    cur_cutoff = maxval(this%a(ti,:)*this%sigma(ti,:))
    n_neigh_i = n_neighbours(at, i)

      do ji=1, n_neigh_i
      j = neighbour(at, i, ji, drij_mag, cosines = drij, max_dist = cur_cutoff)

      if(associated(atom_mask_pointer) .and. atom_mask_exclude_all) then
         if(.not. atom_mask_pointer(j)) cycle
      endif

      if (j <= 0) cycle
      if (drij_mag .feq. 0.0_dp) cycle

      tj = get_type(this%type_of_atomic_num, at%Z(j))

      do ki=1, n_neigh_i

        if (ki <= ji) cycle

        k = neighbour(at, i, ki, drik_mag, cosines = drik, max_dist=cur_cutoff)

        if(associated(atom_mask_pointer) .and. atom_mask_exclude_all) then
           if(.not. atom_mask_pointer(k)) cycle
        endif

        if (k <= 0) cycle
        if (drik_mag .feq. 0.0_dp) cycle

        tk = get_type(this%type_of_atomic_num, at%Z(k))
        if (drik_mag > this%d2(ti, tj, tk)) cycle
        if (drij_mag > this%d1(ti, tj, tk)) cycle

        drij_dot_drik = sum(drij*drik)
        if (present(e) .or. present(local_e)) then
          de = this%eps3(ti,tj,tk)*f3(this, drij_mag, drik_mag, drij_dot_drik, ti, tj, tk)

          if (present(local_e)) then
#ifdef _OPENMP
            private_local_e(i) = private_local_e(i) + de
#else
            local_e(i) = local_e(i) + de
#endif
          endif
          if (present(e)) then
#ifdef _OPENMP
            private_e = private_e + de*w_f
#else
            e = e + de*w_f
#endif
          endif
        endif
        if (present(f) .or. present(virial)) then
          call df3_dr(this, drij_mag, drik_mag, drij_dot_drik, ti, tj, tk, de_drij, de_drik, de_dcos_ijk)
          drij_dri = drij
          drij_drj = -drij !grad.|vec(r_i) -vec(r_j)|= -( vec(r_i) -vec(r_j))/|vec(r_i) -vec(r_j)| = -drij
          drik_dri = drik
          drik_drk = -drik

! dcos_ijk = drij_dot_drik
! (ri_x - rj_x)*(ri_x - rk_x) + (ri_y - rj_y)*(ri_y-rk_y) / (rij rik)
!
! d/dri
!
! ((rij rik)(rij_x + rik_x) - (rij . rik)(rij rik_x / rik + rik rij_x/rij)) / (rij rik)^2
!
! rij rik ( rij_x + rik_x) -  (rij.rik)(rij rik_x / rik + rik rij_x/rij)
! ---------------------------------------------------------------------
!                       rij^2 rik^2
!
! rij_x + rik_x     cos_ijk (rij rik_x / rik + rik rij_x /rij)
! -------------  - ------------------------------------------
! rij rik                    rij rik
!
! rij_x + rik_x
! ------------- - cos_ijk (rik_x / rik^2 + rij_x / rij^2)
!    rij rik
!
!
! rhatij_x/rik + rhatik_x/rij - cos_ijk(rhatik_x/rik + rhatij_x/rij)
!
!
! d/drj
!
! ((rij rik)(-rik_x) - (rij.rik)(rik (-rij_x)/rij)) / (rij rik)^2
!
!  -rik_x     (rij.rik) ( rik rij_x/rij)
! -------  + ---------------------------
! rij rik           rij^2 rik^2
!
! -rhatik_x/rij + cos_ijk rhatij_x/rij
!
! d/drij
!
! (ri_x - rj_x)*(ri_x - rk_x) + (ri_y - rj_y)*(ri_y-rk_y) / (rij rik)
!
! (rij rik) rik_x - (rij . rik) (rij_x/rij rik)
! ---------------------------------------------
!                    (rij rik)^2
!
! rik_x      (rij.rik) rij_x
! ------   - ---------------
! rij rik    rij^3 rik
!
! rhatik_x      (rhatij . rhatik) rhatij_x
! --------   -  ------------------------
! rij           rij
!
! rhatik_x      cos_ijk * rhatij_x
! --------   -  ------------------
!   rij                 rij

          dcos_ijk_dri = drij/drik_mag + drik/drij_mag - drij_dot_drik * (drik/drik_mag + drij/drij_mag)
          dcos_ijk_drj = -drik/drij_mag + drij_dot_drik * drij/drij_mag
          dcos_ijk_drk = -drij/drik_mag + drij_dot_drik * drik/drik_mag

          if (present(f)) then
#ifdef _OPENMP
            private_f(:,i) = private_f(:,i) + w_f*this%eps3(ti,tj,tk)*(de_drij*drij_dri(:) + de_drik*drik_dri(:) + &
                                                       de_dcos_ijk * dcos_ijk_dri(:))
            private_f(:,j) = private_f(:,j) + w_f*this%eps3(ti,tj,tk)*(de_drij*drij_drj(:) + de_dcos_ijk*dcos_ijk_drj(:))
            private_f(:,k) = private_f(:,k) + w_f*this%eps3(ti,tj,tk)*(de_drik*drik_drk(:) + de_dcos_ijk*dcos_ijk_drk(:))
#else
            f(:,i) = f(:,i) + w_f*this%eps3(ti,tj,tk)*(de_drij*drij_dri(:) + de_drik*drik_dri(:) + &
                                                       de_dcos_ijk * dcos_ijk_dri(:))
            f(:,j) = f(:,j) + w_f*this%eps3(ti,tj,tk)*(de_drij*drij_drj(:) + de_dcos_ijk*dcos_ijk_drj(:))
            f(:,k) = f(:,k) + w_f*this%eps3(ti,tj,tk)*(de_drik*drik_drk(:) + de_dcos_ijk*dcos_ijk_drk(:))
#endif
          end if

          if( present(virial) .or. present(local_virial) ) virial_i = &
              w_f*this%eps3(ti,tj,tk)*( &
              de_drij*(drij .outer. drij)*drij_mag + de_drik*(drik .outer. drik)*drik_mag + &
              de_dcos_ijk * ((drik .outer. drij) - drij_dot_drik * (drij .outer. drij)) + &
              de_dcos_ijk * ((drij .outer. drik) - drij_dot_drik * (drik .outer. drik)) )

          if (present(virial)) then
#ifdef _OPENMP
            private_virial = private_virial - virial_i
#else
            virial = virial - virial_i
#endif
          end if
          if (present(local_virial)) then
#ifdef _OPENMP
            private_local_virial(:,i) = private_local_virial(:,i) - reshape(virial_i,(/9/))
#else
            local_virial(:,i) = local_virial(:,i) - reshape(virial_i,(/9/))
#endif
          end if
        endif

      end do

    end do
  end do

#ifdef _OPENMP
!$omp critical
  if (present(e)) e = e + private_e
  if (present(f)) f = f + private_f
  if (present(local_e)) local_e = local_e + private_local_e
  if (present(virial)) virial = virial + private_virial
  if (present(local_virial)) local_virial = local_virial + private_local_virial
!$omp end critical

  if(allocated(private_f)) deallocate(private_f)
  if(allocated(private_local_e)) deallocate(private_local_e)
  if(allocated(private_local_virial)) deallocate(private_local_virial)

!$omp end parallel
#endif

  if (present(mpi)) then
     if (present(e)) e = sum(mpi, e)
     if (present(local_e)) call sum_in_place(mpi, local_e)
     if (present(f)) call sum_in_place(mpi, f)
     if (present(virial)) call sum_in_place(mpi, virial)
     if (present(local_virial)) call sum_in_place(mpi, local_virial)
  endif

end subroutine IPModel_SW_VP_Calc

function f3(this, rij_i, rik_i, cos_ijk, ti, tj, tk)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: rij_i, rik_i, cos_ijk
  integer, intent(in) :: ti, tj, tk
  real(dp) :: f3

  real(dp) :: rij, rik
  !We give all cosines as parameters from xml file

  rij = rij_i
  rik = rik_i

  if (rij_i >= this%d1(ti,tj,tk) .or. rik_i >= this%d2(ti,tj,tk)) then
     f3 = 0.0_dp
    return
  endif

  f3 = this%lambda(ti,tj,tk)*exp( this%gamma1(ti,tj,tk)/(rij-this%d1(ti,tj,tk)) + this%gamma2(ti,tj,tk)/(rik-this%d2(ti,tj,tk)) )* &
     (cos_ijk - this%costheta0(ti,tj,tk))**2

end function f3

subroutine df3_dr(this, rij_i, rik_i, cos_ijk, ti, tj, tk, de_drij, de_drik, de_dcos_ijk)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: rij_i, rik_i, cos_ijk
  integer, intent(in) :: ti, tj, tk
  real(dp), intent(out) :: de_drij, de_drik, de_dcos_ijk

  real(dp) :: rij, rik
  real(dp) :: expf, cosf

  rij = rij_i
  rik = rik_i

  if (rij_i >= this%d1(ti,tj,tk) .or. rik_i >= this%d2(ti,tj,tk)) then ! Anke_new
    de_drij = 0.0_dp
    de_drik = 0.0_dp
    de_dcos_ijk = 0.0_dp
    return
  endif

  expf = this%lambda(ti,tj,tk)*exp(this%gamma1(ti,tj,tk)/(rij - this%d1(ti,tj,tk)) + this%gamma2(ti,tj,tk)/(rik - this%d2(ti,tj,tk)))
  cosf = (cos_ijk -this%costheta0(ti,tj,tk))

  de_drij = -expf * this%gamma1(ti,tj,tk)/(rij-this%d1(ti,tj,tk))**2 * cosf**2
  de_drik = -expf * this%gamma2(ti,tj,tk)/(rik-this%d2(ti,tj,tk))**2 * cosf**2
  de_dcos_ijk = expf * 2.0_dp * cosf

end subroutine df3_dr

!************Begin: Anke modification**********

function f2SiSi(this, ri, ti, tj)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: ri
  integer, intent(in) :: ti, tj

  real(dp) :: f2SiSi

!  if (ri >= this%a(ti,tj)) then
!    f2SiSi = 0.0_dp
!    return
!  endif

!  f2SiSi = this%AA(ti,tj)*(1.0_dp+3.2_dp*qi*qj)*faqj*faqi*(this%BB(ti,tj) * r**(-this%p(ti,tj)) - r**(-this%q(ti,tj))) * exp(1.0_dp/(r-this%a(ti,tj)))
! for now, fixed coordination. If Si is fourfold coordinated, in Si0_2 is faq(Si) = 0.
!
! Potential has to be modified if one wants Si-Si interactions as well

  f2SiSi = 0.0_dp

end function f2SiSi

function df2SiSi_dr(this, ri, ti, tj)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: ri
  integer, intent(in) :: ti, tj
  real(dp) :: df2SiSi_dr

  real(dp) :: expf, dexpf_dr, powf, dpowf_dr

 ! if (r >= this%a(ti,tj)) then
 !   df2SiSi_dr = 0.0_dp

 !   return
 ! endif

 ! expf = exp(this%sigma(ti,tj)/(r-this%a(ti,tj)))
 ! dexpf_dr = -expf*this%sigma(ti,tj)/(r-this%a(ti,tj))**2

 ! powf = this%BB(ti,tj)*(r**(-this%p(ti,tj))) - r**(-this%q(ti,tj))
 ! dpowf_dr = -this%p(ti,tj)*this%BB(ti,tj)*(r**(-this%p(ti,tj)-1)) + this%q(ti,tj)*r**(-this%q(ti,tj)-1)

   df2SiSi_dr = 0.0_dp

end function df2SiSi_dr

function f2SiO(this, r, ti, tj)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: r
  integer, intent(in) :: ti, tj
  real(dp) :: f2SiO

  real(dp) :: expf
  real(dp) :: f2SiO_at_cutoff, df2SiO_dr_at_cutoff

! Cutoff sigma*a for longrange non-coulombic terms

  if (r >= (this%sigma(ti,tj)*this%a(ti,tj))) then
    f2SiO = 0.0_dp
    return
  endif

  expf = exp(-r/this%a(ti,tj))
  f2SiO_at_cutoff = (this%C_0(ti,tj)-this%C_1(ti,tj)*1.6_dp)*(this%sigma(ti,tj)*this%a(ti,tj))**(-9.0_dp)-this%D_SiO(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj))**(-4.0_dp)*exp(-this%sigma(ti,tj))
  df2SiO_dr_at_cutoff= -9.0_dp*(this%C_0(ti, tj)-this%C_1(ti,tj)*1.6_dp)*(this%sigma(ti,tj)*this%a(ti,tj)) **(-10.0_dp) &
                         + 4.0_dp*this%D_SiO(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj)) **(-5.0_dp)* exp(-this%sigma(ti,tj)) &
                         + this%D_SiO(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj) )**(-4.0_dp)*exp(-this%sigma(ti,tj))/this%a(ti,tj)

  f2SiO = (this%C_0(ti,tj)-this%C_1(ti,tj)*1.6_dp)*r**(-9.0_dp) - this%D_SiO(ti,tj)*r**(-4.0_dp)*expf  - f2SiO_at_cutoff - (r - this%sigma(ti,tj)*this%a(ti,tj))*df2SiO_dr_at_cutoff

end function f2SiO

function df2SiO_dr(this, r, ti, tj)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: r
  integer, intent(in) :: ti, tj
  real(dp) :: df2SiO_dr
  real(dp) :: expf
  real(dp) :: df2SiO_dr_at_cutoff

  if (r >= (this%sigma(ti,tj)*this%a(ti,tj))) then
    df2SiO_dr = 0.0_dp
    return
  endif

  expf = exp(-r/this%a(ti,tj))
  df2SiO_dr_at_cutoff= -9.0_dp*(this%C_0(ti, tj)-this%C_1(ti,tj)*1.6_dp)*(this%sigma(ti,tj)*this%a(ti,tj)) **(-10.0_dp) &
                         + 4.0_dp*this%D_SiO(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj)) **(-5.0_dp)* exp(-this%sigma(ti,tj)) &
                         + this%D_SiO(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj) )**(-4.0_dp)*exp(-this%sigma(ti,tj))/this%a(ti,tj)

  df2SiO_dr = -9.0_dp*(this%C_0(ti, tj)-this%C_1(ti,tj)*1.6_dp)*r**(-10.0_dp) + 4.0_dp*this%D_SiO(ti,tj)*r**(-5.0_dp)*expf + this%D_SiO(ti,tj)*r**(-4.0_dp)*expf/this%a(ti,tj) - df2SiO_dr_at_cutoff

end function df2SiO_dr

function f2OO(this, r, ti, tj)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: r
  integer, intent(in) :: ti, tj
  real(dp) :: f2OO
  real(dp) :: expf
  real(dp) :: f2OO_at_cutoff, df2OO_dr_at_cutoff

  if (r >= (this%sigma(ti,tj)*this%a(ti,tj))) then
    f2OO = 0.0_dp
    return
  endif

  expf = exp(-r/this%a(ti,tj))
  f2OO_at_cutoff= (this%C_OOprime(ti, tj))*(this%sigma(ti,tj)*this%a(ti,tj))**(-7.0_dp) - this%D_OOprime(ti, tj)*(this%sigma(ti,tj)*this%a(ti,tj))**(-4.0_dp)*exp(-this%sigma(ti,tj))
 ! Cutoff is at r = sigma*a, thus can be altered via parameter list
  df2OO_dr_at_cutoff=  -7.0_dp*(this%C_OOprime(ti,tj))*(this%sigma(ti,tj)*this%a(ti,tj) )**(-8.0_dp) + 4.0_dp*this%D_OOprime(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj) )**(-5.0_dp)* exp(-this%sigma(ti,tj)) &
                       +  this%D_OOprime(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj) )**(-4.0_dp)*exp(-this%sigma(ti,tj))/this%a(ti,tj)

  f2OO = (this%C_OOprime(ti, tj))*r**(-7.0_dp) - this%D_OOprime(ti, tj)*r**(-4.0_dp)*expf  - f2OO_at_cutoff &
        - (r - this%sigma(ti,tj)*this%a(ti,tj))*df2OO_dr_at_cutoff

end function f2OO

function df2OO_dr(this, r, ti, tj)
  type(IPModel_SW_VP), intent(in) :: this
  real(dp), intent(in) :: r
  integer, intent(in) :: ti, tj
  real(dp) :: df2OO_dr, df2OO_dr_at_cutoff

  real(dp) :: expf

  if (r >= (this%sigma(ti,tj)*this%a(ti,tj))) then
   df2OO_dr = 0.0_dp
    return
  endif

  expf = exp(-r/this%a(ti,tj))

  df2OO_dr_at_cutoff = -7.0_dp*(this%C_OOprime(ti,tj))*(this%sigma(ti,tj)*this%a(ti,tj) )**(-8.0_dp) + 4.0_dp*this%D_OOprime(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj) )**(-5.0_dp)* exp(-this%sigma(ti,tj)) &
                       +  this%D_OOprime(ti,tj)*(this%sigma(ti,tj)*this%a(ti,tj) )**(-4.0_dp)*exp(-this%sigma(ti,tj))/this%a(ti,tj)

  df2OO_dr = -7.0_dp*(this%C_OOprime(ti,tj))*r**(-8.0_dp) + 4.0_dp*this%D_OOprime(ti,tj)*r**(-5.0_dp)*expf +  this%D_OOprime(ti,tj)*r**(-4.0_dp)*expf/this%a(ti,tj)- df2OO_dr_at_cutoff

end function df2OO_dr

!*************End: Anke modification**********

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!% XML param reader functions.
!% An example for the input file.
!%
!%> <SW_VP_potential>
!%>  <Potential label="SW_VP_Potential" init_args="Sum init_args_pot1={IP SW_VP} init_args_pot2={IP Coulomb}"/>
!%>
!%>   <Coulomb_params n_types="3" cutoff="6.0" method="ewald" ewald_error="1.0e-3" label="default_ewald">
!%>   <!-- <per_type_data type="1" atomic_num="1" charge="1.0"/>
!%>      <per_type_data type="2" atomic_num="8" charge="-0.8"/>
!%>      <per_type_data type="3" atomic_num="14" charge="1.6"/> charge can be directly read from extended xyz file -->
!%>   </Coulomb_params>
!%>
!%>
!%>  <SW_VP_params n_types="3" label="PRB_31_plus_H">
!%>  <comment> Danny potential J ChemPhys 127, 204704(2007).The potential has to be called as a pot_init_args='Sum init_args_pot1={IP SW_VP} init_args_pot2={IP Coulomb}'</comment>
!%>  <per_type_data type="1" atomic_num="1" />
!%>  <per_type_data type="2" atomic_num="8" />
!%>  <per_type_data type="3" atomic_num= "14"/>
!%>  <per_pair_data atnum_i="1" atnum_j="1" AA="0.0" BB="0.0"
!%>        p="0" q="0" a="1.0" sigma="1.0" eps="0.0" C_0="0.0" C_1="0.0" D_SiO="0" C_OOprime="0" D_OOprime="0" />
!%>  <per_pair_data atnum_i="1" atnum_j="8" AA="0.0" BB="0.0"
!%>        p="0" q="0" a="1.0" sigma="1.0" eps="0.0" C_0="0.0" C_1="0.0" D_SiO="0" C_OOprime="0" D_OOprime="0" />
!%>  <per_pair_data atnum_i="1" atnum_j="14" AA="8.581214" BB="0.0327827"
!%>        p="4" q="0" a="1.25" sigma="2.537884" eps="0.000" C_0="0.0" C_1="0.0" D_SiO="0" C_OOprime="0" D_OOprime="0"/>
!%>  <per_pair_data atnum_i="8" atnum_j="8" AA="0.0" BB="0.0"
!%>        p="0" q="0" a="4.43" sigma="1.24" eps="14.39" C_0="0.0" C_1="0.0" D_SiO="0" C_OOprime="51.692" D_OOprime="1.536"/>
!%>  <per_pair_data atnum_i="8" atnum_j="14" AA="0.0" BB="0.0"
!%>        p="0" q="0" a="4.43" sigma="1.24" eps="14.39" C_0="14.871" C_1="2.178" D_SiO="3.072" C_OOprime="0.0" D_OOprime="0.0"/>
!%>  <per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
!%>        p="4" q="0" a="1.80" sigma="2.0951" eps="0" C_0="0.0" C_1="0.0" D_SiO="0" C_OOprime="0.0" D_OOprime="0.0"/>
!%>
!%>  <!-- triplet terms: atnum_c is the center atom, neighbours j and k -->
!%>  <per_triplet_data atnum_c="1"  atnum_j="1"  atnum_k="1"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="1"  atnum_j="1"  atnum_k="8"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="1"  atnum_j="1"  atnum_k="14"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="1"  atnum_j="8"  atnum_k="8"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="1"  atnum_j="8"  atnum_k="14"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="1"  atnum_j="14" atnum_k="14"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="8" atnum_j="1" atnum_k="1"
!%>        lambda="3.46" gamma1="0.48" gamma2="0.48" d1="1.24" d2="1.24"  costheta0="-0.5373"  eps="0.0" />
!%>  <per_triplet_data atnum_c="8" atnum_j="1" atnum_k="8"
!%>        lambda="3.46" gamma1="0.48" gamma2="0.48" d1="1.24" d2="1.24"  costheta0="-0.5373"  eps="0.0" />
!%>  <per_triplet_data atnum_c="8" atnum_j="1" atnum_k="14"
!%>        lambda="0.521" gamma1="1" gamma2="1" d1="2.6" d2="2.6"  costheta0="-0.5373"  eps="14.39" />
!%>  <per_triplet_data atnum_c="8" atnum_j="8" atnum_k="8"
!%>        lambda="3.46" gamma1="0.48" gamma2="0.48" d1="1.24" d2="1.24"  costheta0="-0.5373"  eps="0.0" />
!%>  <per_triplet_data atnum_c="8" atnum_j="8" atnum_k="14"
!%>        lambda="3.46" gamma1="0.48" gamma2="0.48" d1="1.24" d2="1.24"  costheta0="-0.5373"  eps="0.0" />
!%>  <per_triplet_data atnum_c="8" atnum_j="14" atnum_k="14"
!%>        lambda="1.4" gamma1="1" gamma2="1" d1="2.6" d2="2.6"  costheta0="-0.777"  eps="14.39" />
!%>  <per_triplet_data atnum_c="14" atnum_j="1"  atnum_k="1"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="14" atnum_j="1"  atnum_k="8"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="14" atnum_j="1"  atnum_k="14"
!%>        lambda="21.0" gamma1="1.20" gamma2="1.20" d1="1.80" d2="1.80" costheta0="-0.333" eps="0.0" />
!%>  <per_triplet_data atnum_c="14" atnum_j="8"  atnum_k="8"
!%>        lambda="0.350" gamma1="1" gamma2="1" d1="2.6" d2="2.6" costheta0="-0.333" eps="14.39" />
!%>  <!--In original paper, the c=14,j=14,k=8 term is listed, hence d2 and d1 and gamma1 and gamma2 are interchanged  -->
!%>  <per_triplet_data atnum_c="14" atnum_j="8" atnum_k="14"
!%>        lambda="3.164" gamma1="4.06" gamma2="0.52" d1="3.981" d2="2.933"  costheta0="-0.333" eps="14.39" />
!%>  <per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
!%>        lambda="3.164" gamma1="2.51" gamma2="2.51" d1="3.771" d2="3.771" costheta0="-0.333"  eps="14.39" />
!%>  </SW_VP_params>
!%>
!%>  <constraints N="1">
!%> <!-- All OH-bond lengths have to be constrained to avoid coulomb clash between OH for surfaces-->
!%>  </constraints>
!%> </SW_VP_potential>
!%
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

subroutine IPModel_startElement_handler(URI, localname, name, attributes)
  character(len=*), intent(in)   :: URI
  character(len=*), intent(in)   :: localname
  character(len=*), intent(in)   :: name
  type(dictionary_t), intent(in) :: attributes

  integer status
  character(len=1024) :: value
  integer ti, tj, tk, Zi, Zj, Zk

  if (name == 'SW_VP_params') then ! new SW_VP stanza
    call print("startElement_handler SW_VP_params", PRINT_NERD)
    if (parse_in_ip) &
      call system_abort("IPModel_startElement_handler entered SW_VP_params with parse_in true. Probably a bug in FoX (4.0.1, e.g.)")

    if (parse_matched_label) then
       call print("SW_VP_params startElement_handler bailing because we already matched our label", PRINT_NERD)
       return ! we already found an exact match for this label
    end if

    call QUIP_FoX_get_value(attributes, 'label', value, status)
    if (status /= 0) value = ''

    call print("SW_VP_params startElement_handler found xml label '"//trim(value)//"'", PRINT_NERD)

    if (len(trim(parse_ip%label)) > 0) then ! we were passed in a label
      call print("SW_VP_params startElement_handler was passed in label '"//trim(parse_ip%label)//"'", PRINT_NERD)
      if (value == parse_ip%label) then ! exact match
        parse_matched_label = .true.
        parse_in_ip = .true.
      else ! no match
        call print("SW_VP_params startElement_handler got label didn't match", PRINT_NERD)
        parse_in_ip = .false.
      endif
    else ! no label passed in
      call print("SW_VP_params startElement_handler was not passed in a label", PRINT_NERD)
      parse_in_ip = .true.
    endif

    if (parse_in_ip) then
      if (parse_ip%n_types /= 0) then
        call print("SW_VP_params startElement_handler finalising old data, restarting to parse new section", PRINT_NERD)
        call finalise(parse_ip)
      endif

      call QUIP_FoX_get_value(attributes, "n_types", value, status)
      if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find n_types")
      read (value, *) parse_ip%n_types

      allocate(parse_ip%atomic_num(parse_ip%n_types))
      parse_ip%atomic_num = 0
      allocate(parse_ip%a(parse_ip%n_types,parse_ip%n_types))
      parse_ip%a = 0.0_dp
      allocate(parse_ip%AA(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%BB(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%p(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%q(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%sigma(parse_ip%n_types,parse_ip%n_types))
      parse_ip%sigma = 0.0_dp
      allocate(parse_ip%eps2(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%C_0(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%C_1(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%D_SiO(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%C_OOprime(parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%D_OOprime(parse_ip%n_types,parse_ip%n_types))

      allocate(parse_ip%lambda(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%gamma1(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%gamma2(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%eps3(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%d1(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%d2(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
      allocate(parse_ip%costheta0(parse_ip%n_types,parse_ip%n_types,parse_ip%n_types))
    endif

  elseif (parse_in_ip .and. name == 'per_type_data') then

    call QUIP_FoX_get_value(attributes, "type", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find type")
    read (value, *) ti

    call QUIP_FoX_get_value(attributes, "atomic_num", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find atomic_num")
    read (value, *) parse_ip%atomic_num(ti)

    if (allocated(parse_ip%type_of_atomic_num)) deallocate(parse_ip%type_of_atomic_num)
    allocate(parse_ip%type_of_atomic_num(maxval(parse_ip%atomic_num)))
    parse_ip%type_of_atomic_num = 0
    do ti=1, parse_ip%n_types
      if (parse_ip%atomic_num(ti) > 0) &
        parse_ip%type_of_atomic_num(parse_ip%atomic_num(ti)) = ti
    end do

  elseif (parse_in_ip .and. name == 'per_pair_data') then

    call QUIP_FoX_get_value(attributes, "atnum_i", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find atnum_i")
    read (value, *) Zi
    call QUIP_FoX_get_value(attributes, "atnum_j", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find atnum_j")
    read (value, *) Zj

    ti = get_type(parse_ip%type_of_atomic_num,Zi)
    tj = get_type(parse_ip%type_of_atomic_num,Zj)

    call QUIP_FoX_get_value(attributes, "AA", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find AA")
    read (value, *) parse_ip%AA(ti,tj)
    call QUIP_FoX_get_value(attributes, "BB", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find BB")
    read (value, *) parse_ip%BB(ti,tj)
    call QUIP_FoX_get_value(attributes, "p", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find p")
    read (value, *) parse_ip%p(ti,tj)
    call QUIP_FoX_get_value(attributes, "q", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find q")
    read (value, *) parse_ip%q(ti,tj)
    call QUIP_FoX_get_value(attributes, "a", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find a")
    read (value, *) parse_ip%a(ti,tj)
    call QUIP_FoX_get_value(attributes, "sigma", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find sigma")
    read (value, *) parse_ip%sigma(ti,tj)
    call QUIP_FoX_get_value(attributes, "eps", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find eps")
    read (value, *) parse_ip%eps2(ti,tj)
    call QUIP_FoX_get_value(attributes, "C_0", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find C_0")
    read (value, *) parse_ip%C_0(ti,tj)
    call QUIP_FoX_get_value(attributes, "C_1", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find C_1")
    read (value, *) parse_ip%C_1(ti,tj)
    call QUIP_FoX_get_value(attributes, "D_SiO", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find D_SiO")
    read (value, *) parse_ip%D_SiO(ti,tj)
    call QUIP_FoX_get_value(attributes, "C_OOprime", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find C_OOprime")
    read (value, *) parse_ip%C_OOprime(ti,tj)
    call QUIP_FoX_get_value(attributes, "D_OOprime", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find D_OOprime")
    read (value, *) parse_ip%D_OOprime(ti,tj)

    if (ti /= tj) then
      parse_ip%AA(tj,ti) = parse_ip%AA(ti,tj)
      parse_ip%BB(tj,ti) = parse_ip%BB(ti,tj)
      parse_ip%p(tj,ti) = parse_ip%p(ti,tj)
      parse_ip%q(tj,ti) = parse_ip%q(ti,tj)
      parse_ip%a(tj,ti) = parse_ip%a(ti,tj)
      parse_ip%sigma(tj,ti) = parse_ip%sigma(ti,tj)
      parse_ip%eps2(tj,ti) = parse_ip%eps2(ti,tj)
      parse_ip%C_0(tj,ti) = parse_ip%C_0(ti,tj)
      parse_ip%C_1(tj,ti) = parse_ip%C_1(ti,tj)
      parse_ip%D_SiO(tj,ti) = parse_ip%D_SiO(ti,tj)
      parse_ip%C_OOprime(tj,ti) = parse_ip%C_OOprime(ti,tj)
      parse_ip%D_OOprime(tj,ti) = parse_ip%D_OOprime(ti,tj)
    endif

    parse_ip%cutoff = maxval(parse_ip%a*parse_ip%sigma)

  elseif (parse_in_ip .and. name == 'per_triplet_data') then

    call QUIP_FoX_get_value(attributes, "atnum_c", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find atnum_c")
    read (value, *) Zi
    call QUIP_FoX_get_value(attributes, "atnum_j", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find atnum_j")
    read (value, *) Zj
    call QUIP_FoX_get_value(attributes, "atnum_k", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find atnum_k")
    read (value, *) Zk

    ti = get_type(parse_ip%type_of_atomic_num,Zi)
    tj = get_type(parse_ip%type_of_atomic_num,Zj)
    tk = get_type(parse_ip%type_of_atomic_num,Zk)

    call QUIP_FoX_get_value(attributes, "lambda", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find lambda")
    read (value, *) parse_ip%lambda(ti,tj,tk)
    call QUIP_FoX_get_value(attributes, "gamma1", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find gamma1")
    read (value, *) parse_ip%gamma1(ti,tj,tk)
    call QUIP_FoX_get_value(attributes, "gamma2", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find gamma2")
    read (value, *) parse_ip%gamma2(ti,tj,tk)
    call QUIP_FoX_get_value(attributes, "d1", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find d1")
    read (value, *) parse_ip%d1(ti,tj,tk)
    call QUIP_FoX_get_value(attributes, "d2", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find d2")
    read (value, *) parse_ip%d2(ti,tj,tk)
    call QUIP_FoX_get_value(attributes, "costheta0", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find costheta0")
    read (value, *) parse_ip%costheta0(ti,tj,tk)
    call QUIP_FoX_get_value(attributes, "eps", value, status)
    if (status /= 0) call system_abort ("IPModel_SW_VP_read_params_xml cannot find eps")
    read (value, *) parse_ip%eps3(ti,tj,tk)

    if (tj /= tk) then
      parse_ip%lambda(ti,tk,tj) = parse_ip%lambda(ti,tj,tk)
      parse_ip%gamma1(ti,tk,tj) = parse_ip%gamma2(ti,tj,tk)     ! ti is central atom for tj /=tk, the legs change and hence gamma1 and gamma2 interchange
      parse_ip%gamma2(ti,tk,tj) = parse_ip%gamma1(ti,tj,tk)     !
      parse_ip%d1(ti,tk,tj) = parse_ip%d2(ti,tj,tk)             !  ti is central atom for tj /=tk, the legs change and hence d1 and d2 interchange
      parse_ip%d2(ti,tk,tj) = parse_ip%d1(ti,tj,tk)             !
      parse_ip%eps3(ti,tk,tj) = parse_ip%eps3(ti,tj,tk)
      parse_ip%costheta0(ti,tk,tj) = parse_ip%costheta0(ti,tj,tk)
    endif
  endif

end subroutine IPModel_startElement_handler

subroutine IPModel_endElement_handler(URI, localname, name)
  character(len=*), intent(in)   :: URI
  character(len=*), intent(in)   :: localname
  character(len=*), intent(in)   :: name

  if (parse_in_ip) then
    if (name == 'SW_VP_params') then
      call print("endElement_handler SW_VP_params", PRINT_NERD)
      parse_in_ip = .false.
    endif
  endif

end subroutine IPModel_endElement_handler

subroutine IPModel_SW_VP_read_params_xml(this, param_str)
  type(IPModel_SW_VP), intent(inout), target :: this
  character(len=*) :: param_str

  type (xml_t) :: fxml

  if (len(trim(param_str)) <= 0) return

  parse_ip => this
  parse_in_ip = .false.
  parse_matched_label = .false.

  call open_xml_string(fxml, param_str)

  call parse(fxml, &
    startElement_handler = IPModel_startElement_handler, &
    endElement_handler = IPModel_endElement_handler)

  call close_xml_t(fxml)

  if (this%n_types == 0) then
    call system_abort("IPModel_SW_VP_read_params_xml parsed file, but n_types = 0")
  endif

    end subroutine IPModel_SW_VP_read_params_xml

!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
!X
!X printing
!% Printing of SW_VP parameters: number of different types, cutoff radius, atomic numbers, ect.
!X
!XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

subroutine IPModel_SW_VP_Print (this, file)
  type(IPModel_SW_VP), intent(in) :: this
  type(Inoutput), intent(inout),optional :: file

  integer ti, tj, tk

  call Print("IPModel_SW_VP : Combined Stillinger-Weber/Vashishta potential", file=file)
  call Print("IPModel_SW_VP : n_types = " // this%n_types // " cutoff = " // this%cutoff, file=file)

  do ti=1, this%n_types
    call Print ("IPModel_SW_VP : type "// ti // " atomic_num " // this%atomic_num(ti), file=file)
    call verbosity_push_decrement()
    do tj=1, this%n_types
      call Print ("IPModel_SW_VP : pair interaction ti tj " // ti // " " // tj // " Zi Zj " // this%atomic_num(ti) // &
        " " // this%atomic_num(tj), file=file)
      call Print ("IPModel_SW_VP : pair " // this%AA(ti,tj) // " " // this%BB(ti,tj) // " " // this%p(ti,tj) // " " // &
        this%q(ti,tj) // " " // this%a(ti,tj), file=file)
      call Print ("IPModel_SW_VP :      " // this%sigma(ti,tj) // " " // this%eps2(ti,tj), file=file)
      do tk=1, this%n_types
        call Print ("IPModel_SW_VP : triplet interaction ti tj " // ti // " " // tj // " " // tk // &
          " Zi Zj Zk " // this%atomic_num(ti) // " " // this%atomic_num(tj) // " " // this%atomic_num(tk), file=file)
        call Print ("IPModel_SW_VP : triplet " // this%lambda(ti,tj,tk) // " " // this%gamma1(ti,tj,tk) // " " // &
          this%eps3(ti,tj,tk), file=file)
      end do
    end do
    call verbosity_pop()
  end do

end subroutine IPModel_SW_VP_Print

end module IPModel_SW_VP_module

max_stars_repo_name: Sideboard/QUIP

max_issues_repo_name: Sideboard/QUIP

max_forks_repo_name: Sideboard/QUIP


text:
      subroutine ttime(ip)
c-----------------------------------------------------------------------
c       this subroutine performs the major calculation in this program.
c  it generates a travel time curve corresponding to the present source
c  depth.  the section of the total velocity model of length nwindow
c  formed by ttable is used in an accumulative fashion to calculate the
c  travel times of the upward rays.  downward rays are then calculated
c  by the simple difference from the surface focus travel time curve.
c  branches are flagged as forward, revers, shadow,or upward.  arguments
c  passed are-
c         ip-index to enter surface focus table.
c     inputs coming through common-
c         nwindow-length of window in velocity model for this increment
c         vwork-velocity array giving velocity model window
c         dwork-corresponding depths
c-----------------------------------------------------------------------
	include 'table.common'
	include 'upward.common'
	include 'surface.common'
	include 'control.common'
	include 'window.common'
      character bnorm,revs,shad,up
      parameter(bnorm='d',revs='r',shad='s',up='u')
      parameter(flaglvz=1.0e38)
c--loop on rays.  find the accumulated travel time and distance to this
c--depth.
          do 100 iup=ip+1,np
          rayp=pup(iup)
          call txint(rayp,vwork,dwork,nwindow,tt,xx)
          tup(iup)=tup(iup)+tt
  100     xup(iup)=xup(iup)+xx
c--find smallest p needed
      kount=ip+1
  125 if(xup(kount).lt.xmin) go to 150
      kount=kount + 1
      go to 125
  150 continue
      lastp=kount
c--do not do straight up case twice if kount=np
      if(kount.eq.np) lastp=lastp-1
c--calculate largest allowed ray index
c--pmax=1/vsource
      ipmax=ip+1
c--calculate the reversed segment of upward rays.
c--start with the straight up ray.
      ttab(1)=tup(np)
      xtab(1)=0.
      ptab(1)=0.
      branch(1)=up
      irange = lastp-ipmax+2
c--when low velocity zones are present the pup array is not a
c--continuously decreasing sequence.  when this occurs in the following
c--loop, those points are eliminated by testing the t array.
      itt=2
          do 175 iup=2,irange
          iuprev=lastp-iup+2
c--do not make test for small p.
          if(iuprev.gt.ns) go to 170
               if(t(iuprev).gt.flaglvz) go to 175
  170     ttab(itt)=tup(iuprev)
          xtab(itt)=xup(iuprev)
          ptab(itt)=pup(iuprev)
          branch(itt)=up
          itt=itt+1
  175     continue
      itt=itt-1
c--for sources outside low velocity zones shoot one extra ray
c--straight out.
c--calculate travel times for rays starting out downward
c--first check to see if there are any.  if not exit gracefully.
      if(ipmax.ge.ns) then
          ntab=irange
          return
      endif
      if(t(ip+1).lt.flaglvz) then
          itt=itt+1
          vsource=vwork(nwindow)
          pmax=1./vsource
c--be careful of coasely traced layered models for near surface
c--sources.
          if(ip.le.1) then
               tt = 0.0
               xx = 0.0
          else
               call txint(pmax,v,d,ip,tt,xx)
          endif
c--be careful if the bottom of the window falls on a grid point.
c--(see ttable.)
            if(dwork(nwindow).eq.d(ip)) then
                 dtt=0.0
                 dxx=0.0
            else
                 call txint(pmax,vwork(nwindow-1),dwork(nwindow-1),2,
     $                       dtt,dxx)
            endif
          ttab(itt)=tt+dtt
          xtab(itt)=xx+dxx
          ptab(itt)=pmax
          if(ttab(itt).lt.flaglvz) then
               branch(itt)=up
          else
               branch(itt)=shad
          endif
      endif
          do 200 idown=ipmax,ns
          itt=itt+1
          if(t(idown).gt.flaglvz) then
c--lvz branch
               ttab(itt)=t(idown)
               xtab(itt)=x(idown)
               ptab(itt)=p(idown)
               branch(itt)=shad
          else
c--normal branch
               ttab(itt)=t(idown)-tup(idown)
               xtab(itt)=x(idown)-xup(idown)
               ptab(itt)=p(idown)
c--branch is always named by looking backward one point.
c--first point after a shadow zone is labeled as = shad
               if(xtab(itt-1).gt.flaglvz) then
                    branch(itt) = shad
               elseif(xtab(itt).ge.xtab(itt-1))then
                    branch(itt)=bnorm
               else
                    branch(itt)=revs
               endif
          endif
  200     continue
      ntab=itt
      return
      end

c $Id$

ext: f


text:
great_rsd(w1,l1).
great_rsd(p1,h1).
great_rsd(p1,kk1).
great_rsd(w1,u1).
great_rsd(u1,t1).
great_rsd(aa1,m1).
great_rsd(aa1,x1).
great_rsd(cc1,z1).
great_rsd(a1,kk1).
great_rsd(q1,x1).
great_rsd(e1,u1).
great_rsd(u1,l1).
great_rsd(aa1,ee1).
great_rsd(ee1,h1).
great_rsd(ff1,hh1).
great_rsd(m1,hh1).
great_rsd(w1,ii1).
great_rsd(p1,q1).
great_rsd(k1,dd1).
great_rsd(u1,hh1).
great_rsd(f1,c1).
great_rsd(k1,ii1).
great_rsd(w1,kk1).
great_rsd(kk1,c1).
great_rsd(a1,g1).
great_rsd(t1,h1).
great_rsd(f1,m1).
great_rsd(v1,ii1).
great_rsd(g1,z1).
great_rsd(t1,ee1).
great_rsd(k1,l1).
great_rsd(k1,c1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text:
!   gfortran -cpp -O3 -flto ex_3.f90 -o ex_3
!   ./ex_3

!   "Convergence" testing: remap a smooth (Gaussian) profile
!   between perturbed grids and compute error in the L2-nrm.
!   Re-run with N = 50, 100, 200, ... etc (and N-tmp at 80%)
!   to assess the order of spatial accuracy. With NUll/WENO-
!   LIMIT, all methods should achieve their best-case O(h^p)
!   scaling.
!

#   include "../src/ppr_1d.f90"

    program ex

        use ppr_1d

        implicit none

        integer, parameter :: npos = 50 ! no. edge (old grid)
        integer, parameter :: ntmp = 40 ! no. edge (new grid)
        integer, parameter :: nvar = 1  ! no. variables to remap
        integer, parameter :: ndof = 1  ! no. FV DoF per cell
        integer :: ipos

    !------------------------------ position of cell edges !
        real*8  :: xpos(npos),xtmp(ntmp)
        real*8  :: xdel,xmid

    !-------------------------------- finite-volume arrays !

    !   Arrays represent a "block" of finite-volume tracers
    !   to remap. The 1st dim. is the no. of DoF per cell,
    !   NDOF=1 is a standard finite-volume scheme where the
    !   data is specified as cell means. NDOF>1 is reserved
    !   for future use with DG-style schemes. NVAR is the
    !   number of tracers to remap. Processing tracers in a
    !   batch is typically more efficient than one-by-one.
    !   The last dim. is the no. cells (layers) in the grid.

        real*8  :: init(ndof,nvar,npos-1)
        real*8  :: ftmp(ndof,nvar,ntmp-1)
        real*8  :: fdat(ndof,nvar,npos-1)

    !------------------------------ method data-structures !
        type(rmap_work) :: work
        type(rmap_opts) :: opts
        type(rcon_ends) :: bc_l(nvar)
        type(rcon_ends) :: bc_r(nvar)

    !------------------------------ define a simple domain !

        call linspace(0.d0,1.d0,npos,xpos)
        call linspace(0.d0,1.d0,ntmp,xtmp)

    !------------------------------ setup some simple data !

        do ipos = +1, npos-1

            xmid = xpos(ipos+0) * 0.5d+0 &
    &            + xpos(ipos+1) * 0.5d+0

            init(1,1,ipos) = &
    &   .8d+0 * exp( -75.0d+0 * (xmid - 0.275d+0) ** 2 ) &
    & + .9d+0 * exp(-100.0d+0 * (xmid - 0.500d+0) ** 2 ) &
    & + 1.d+0 * exp(-125.0d+0 * (xmid - 0.725d+0) ** 2 )

        end do

    !------------------------------ specify method options !

        opts%edge_meth = p5e_method     ! 3rd-order edge interp.
        opts%cell_meth = pqm_method     ! PPM method in cells
        opts%cell_lims = null_limit     ! monotone limiter

    !------------------------------ set BC.'s at endpoints !

        bc_l%bcopt = bcon_loose         ! "loose" = extrapolate
        bc_r%bcopt = bcon_loose

    !------------------------------ init. method workspace !

        call work%init(npos,nvar,opts)

    !------------------------------ re-map back-and-forth: !

        fdat = init

        do ipos = +1, +1000

    !------------------------------ re-map from dat-to-tmp !

        call rmap1d(npos,ntmp,nvar,ndof, &
        &           xpos,xtmp,fdat,ftmp, &
        &           bc_l,bc_r,work,opts)

    !------------------------------ re-map from tmp-to-dat !

        call rmap1d(ntmp,npos,nvar,ndof, &
        &           xtmp,xpos,ftmp,fdat, &
        &           bc_l,bc_r,work,opts)

        end do

    !------------------------------ clear method workspace !

        call work%free()

    !------------------------------ dump results to stdout !

        print*,"Cell data: [INIT] [RMAP] "

        do ipos = +1, npos-1

            print *, init(1,1,ipos) &
        &          , fdat(1,1,ipos)

        end do

        print*,"Conservation defect := " &
        &     , sum(init) - sum(fdat)

    !------------------------------ calc. errors in L2 nrm !

        xdel = (xpos(npos)-xpos(1))/(npos-1)

        print*,"Discrete (L2) Error := ",&
        sqrt(sum(xdel * (init(1,1,:) - fdat(1,1,:)) ** 2))

    end program

max_stars_repo_name: mark-petersen/PPR

max_issues_repo_name: mark-petersen/PPR

max_forks_repo_name: mark-petersen/PPR


text:
      SUBROUTINE COLSOL (A,V,MAXA,NN,NWK,NNM,KKK,IOUT)                  COL00001
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . COL00002
C .                                                                   . COL00003
C .   P R O G R A M                                                   . COL00004
C .        TO SOLVE FINITE ELEMENT STATIC EQUILIBRIUM EQUATIONS IN    . COL00005
C .        CORE, USING COMPACTED STORAGE AND COLUMN REDUCTION SCHEME  . COL00006
C .                                                                   . COL00007
C .  - - INPUT VARIABLES - -                                          . COL00008
C .        A(NWK)    = STIFFNESS MATRIX STORED IN COMPACTED FORM      . COL00009
C .        V(NN)     = RIGHT-HAND-SIDE LOAD VECTOR                    . COL00010
C .        MAXA(NNM) = VECTOR CONTAINING ADDRESSES OF DIAGONAL        . COL00011
C .                    ELEMENTS OF STIFFNESS MATRIX IN A              . COL00012
C .        NN        = NUMBER OF EQUATIONS                            . COL00013
C .        NWK       = NUMBER OF ELEMENTS BELOW SKYLINE OF MATRIX     . COL00014
C .        NNM       = NN + 1                                         . COL00015
C .        KKK       = INPUT FLAG                                     . COL00016
C .            EQ. 1   TRIANGULARIZATION OF STIFFNESS MATRIX          . COL00017
C .            EQ. 2   REDUCTION AND BACK-SUBSTITUTION OF LOAD VECTOR . COL00018
C .        IOUT      = UNIT NUMBER USED FOR OUTPUT                    . COL00019
C .                                                                   . COL00020
C .  - - OUTPUT - -                                                   . COL00021
C .        A(NWK)    = D AND L - FACTORS OF STIFFNESS MATRIX          . COL00022
C .        V(NN)     = DISPLACEMENT VECTOR                            . COL00023
C .                                                                   . COL00024
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . COL00025
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)                               COL00026
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . COL00027
C .   THIS PROGRAM IS USED IN SINGLE PRECISION ARITHMETIC ON CRAY     . COL00028
C .   EQUIPMENT AND DOUBLE PRECISION ARITHMETIC ON IBM MACHINES,      . COL00029
C .   ENGINEERING WORKSTATIONS AND PCS. DEACTIVATE ABOVE LINE FOR     . COL00030
C .   SINGLE PRECISION ARITHMETIC.                                    . COL00031
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . COL00032
      DIMENSION A(NWK),V(NN),MAXA(NNM)                                  COL00033
C                                                                       COL00034
C     PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX                COL00035
C                                                                       COL00036
      IF (KKK-2) 40,150,150                                             COL00037
   40 DO 140 N=1,NN                                                     COL00038
      KN=MAXA(N)                                                        COL00039
      KL=KN + 1                                                         COL00040
      KU=MAXA(N+1) - 1                                                  COL00041
      KH=KU - KL                                                        COL00042
      IF (KH) 110,90,50                                                 COL00043
   50 K=N - KH                                                          COL00044
      IC=0                                                              COL00045
      KLT=KU                                                            COL00046
      DO 80 J=1,KH                                                      COL00047
      IC=IC + 1                                                         COL00048
      KLT=KLT - 1                                                       COL00049
      KI=MAXA(K)                                                        COL00050
      ND=MAXA(K+1) - KI - 1                                             COL00051
      IF (ND) 80,80,60                                                  COL00052
   60 KK=MIN0(IC,ND)                                                    COL00053
      C=0.                                                              COL00054
      DO 70 L=1,KK                                                      COL00055
   70 C=C + A(KI+L)*A(KLT+L)                                            COL00056
      A(KLT)=A(KLT) - C                                                 COL00057
   80 K=K + 1                                                           COL00058
   90 K=N                                                               COL00059
      B=0.                                                              COL00060
      DO 100 KK=KL,KU                                                   COL00061
      K=K - 1                                                           COL00062
      KI=MAXA(K)                                                        COL00063
      C=A(KK)/A(KI)                                                     COL00064
      B=B + C*A(KK)                                                     COL00065
  100 A(KK)=C                                                           COL00066
      A(KN)=A(KN) - B                                                   COL00067
  110 IF (A(KN)) 120,120,140                                            COL00068
  120 WRITE (IOUT,2000) N,A(KN)                                         COL00069
      GO TO 800                                                         COL00070
  140 CONTINUE                                                          COL00071
      GO TO 900                                                         COL00072
C                                                                       COL00073
C     REDUCE RIGHT-HAND-SIDE LOAD VECTOR                                COL00074
C                                                                       COL00075
  150 DO 180 N=1,NN                                                     COL00076
      KL=MAXA(N) + 1                                                    COL00077
      KU=MAXA(N+1) - 1                                                  COL00078
      IF (KU-KL) 180,160,160                                            COL00079
  160 K=N                                                               COL00080
      C=0.                                                              COL00081
      DO 170 KK=KL,KU                                                   COL00082
      K=K - 1                                                           COL00083
  170 C=C + A(KK)*V(K)                                                  COL00084
      V(N)=V(N) - C                                                     COL00085
  180 CONTINUE                                                          COL00086
C                                                                       COL00087
C     BACK-SUBSTITUTE                                                   COL00088
C                                                                       COL00089
      DO 200 N=1,NN                                                     COL00090
      K=MAXA(N)                                                         COL00091
  200 V(N)=V(N)/A(K)                                                    COL00092
      IF (NN.EQ.1) GO TO 900                                            COL00093
      N=NN                                                              COL00094
      DO 230 L=2,NN                                                     COL00095
      KL=MAXA(N) + 1                                                    COL00096
      KU=MAXA(N+1) - 1                                                  COL00097
      IF (KU-KL) 230,210,210                                            COL00098
  210 K=N                                                               COL00099
      DO 220 KK=KL,KU                                                   COL00100
      K=K - 1                                                           COL00101
  220 V(K)=V(K) - A(KK)*V(N)                                            COL00102
  230 N=N - 1                                                           COL00103
      GO TO 900                                                         COL00104
C                                                                       COL00105
  800 STOP                                                              COL00106
  900 RETURN                                                            COL00107
C                                                                       COL00108
 2000 FORMAT (//' STOP - STIFFNESS MATRIX NOT POSITIVE DEFINITE',//,    COL00109
     1          ' NONPOSITIVE PIVOT FOR EQUATION ',I8,//,               COL00110
     2          ' PIVOT = ',E20.12 )                                    COL00111
      END                                                               COL00112

ext: for


text:
*...................................................................

      integer function ODS_NVal ( ndim, count )

      implicit NONE

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!
! !ROUTINE:  ODS_NVal
!
! !DESCRIPTION:
!     This routine determines the total number of values written to
!     or read from a NetCDF file.  The number of values is a
!     product of the number for each NetCDF dimension for a variable.
!
! !INTERFACE: nval = ODS_NVal ( ndim, count )
!
! !INPUT PARAMETERS:
      integer ndim         ! Number of NetCDF dimensions
      integer count ( * )  ! Number of values for each dimension
!
! !REVISION HISTORY:
!     13Oct95   C. Redder   Origional version
!
!-------------------------------------------------------------------------

*     Other variables
*     ---------------
      integer nval   ! temporary storage for the number of values
      integer idim   ! index variable for do loop

*     Initialize nval if ndim is reasonable
*     -------------------------------------
      if ( ndim .ge. 1 ) then
         nval = count ( 1 )

*     Otherwise, set numbers of values to zero and return
*     ---------------------------------------------------
      else
         ODS_NVal = 0
         return

      end if

*     Determine nval
*     --------------
      do 10, idim = 2, ndim
         nval = nval * count ( idim )
 10   continue

      ODS_NVal = nval

      return
      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
MODULE genericParams
    !====================================================================
    !Define variables that are jointly used by the global optimization
    !as well as in the objective function.
    !Note, parameters that are specific objective function .
    !should not be here. For them use a separate module.
    !====================================================================
    USE nrtype
    USE stateControl
    IMPLICIT NONE

    INTEGER(I4B) :: p_default_alg = 0 !default algorithm.
          !  0-bobyqa (DFNLS), 1-amoeba (simplex), 2-dfpmin (BFGS Quasi-Newton method)
    INTEGER(I4B) :: p_nx, p_nmom
    INTEGER(I4B) :: p_iprint=2, p_maxeval, p_ninterppt ! BOBYQ variables
    INTEGER(I4B):: p_qr_ndraw, p_maxpoints,p_legitimate
    INTEGER(I4B) :: p_searchType, p_lotteryPoints
    REAL(DP), DIMENSION(:), ALLOCATABLE :: p_wspace
    REAL(DP), DIMENSION(:), ALLOCATABLE :: p_init
    REAL(DP), DIMENSION(:,:), ALLOCATABLE :: p_range, p_bound
    REAL(DP) :: p_fvalmax
    INTEGER(I4B), PARAMETER  :: p_nsmplx=100
        ! initialize the seed for random numbers
    INTEGER(I4B), PARAMETER :: p_SEED=314159265
    REAL(DP), parameter :: p_tolf_amoeba= 1.0d-4    ! Amoeba variables
    INTEGER(I4B), parameter :: p_itmax_amoeba=1000  ! Amoeba variables
    REAL(DP), PARAMETER :: gtol_DFPMIN = 1.0d-06			      ! DFPMIN PARAMETER
    INTEGER,  PARAMETER :: itmax_DFPMIN = 1000      ! DFPMIN PARAMETER
    REAL(DP), parameter :: neg_miss=-10000000.0_dp
    REAL(DP), parameter :: pos_miss= 10000000.0_dp

CONTAINS
    SUBROUTINE initialize(option,seqNo,config)
        !Initialize the instance. This includes parsing the configuration file,
        !setting parameter values, etc.
! 0ptions:
!       -1 = exit state - Stop all instances.
!        0 = cold start - The first invocation of the program, that will set up all the parallel helper files.
!                    Should always be the parameter when this is the first attempt at solving the problem.
!        1 = warm start - after one process is already running, all helper programs should be invoked using
!                    a warm start.
!        2 = update number of Sobol points - update the sobol point parameters over which to search
!                    as well as the local optimizations from these sobol points, but assume everything
!                    else in the config file has not been changed.
!        3 = update number of local minimizations - Increase the number of local minimizations but keep
!                    everything else same in the config file. This option uses the results from previously
!                    found local minimums.
!        4 = Diagnostic - Just evaluate the objective functions once for given initial guess in the config file with
!                    diagnostic option.
!        5 = Local Minimization - Run local minimization once for given initial guess in the config file.

        INTEGER(I4B), INTENT(in) :: option
        INTEGER(I4B), INTENT(out) :: seqNo
        CHARACTER(LEN=25), INTENT(IN), OPTIONAL :: config

        LOGICAL :: THERE_EXISTS,isUpdate
        INTEGER(I4B) :: i, n, fileDesc, openStat
        INTEGER, DIMENSION(:), ALLOCATABLE :: seed

        REAL(DP), DIMENSION(:,:), ALLOCATABLE :: Old_sobolFnVal, Old_SearchResults
        INTEGER(I4B) :: numsolved,old_p_qr_ndraw,legitSobol,old_p_maxpoints
        CHARACTER(LEN=1000) :: errorString

        call random_seed(size = n)
        allocate(seed(n))
        seed(1)=123456
        call random_seed(put = seed)
        deallocate(seed)

        IF(OPTION==4 .or. OPTION==5) THEN
              !not updating, just parse the config file and run the objective value
              !for one set of parameters (option 4).
              !or run the given local minimization around the given point (option 5).
              isUpdate= .FALSE.
              call parseConfig(isUpdate, config)
              close(41)
              open(UNIT=41, FILE='seqNo.dat', STATUS='replace'); write(41,*) 0 ; close(41)
              IF(OPTION==5) THEN
                INQUIRE(FILE='searchResults.dat', EXIST=THERE_EXISTS)
                IF(THERE_EXISTS .eqv. .FALSE.) THEN
                  open(UNIT=41, FILE='searchResults.dat', STATUS='replace'); close(41)
                ENDIF

                INQUIRE(FILE='searchStart.dat', EXIST=THERE_EXISTS)
                IF(THERE_EXISTS .eqv. .FALSE.) THEN
                  open(UNIT=41, FILE='searchStart.dat', STATUS='replace'); close(41)
                ENDIF
              ENDIF
              RETURN
        END IF

        !if config is present,
        IF (PRESENT(config)) THEN

            IF (OPTION == 0) THEN
              ! if not updating, then reset all files
                  !delete any of the files managing the state
              isUpdate= .FALSE.
              close(41)
              open(UNIT=41, FILE='seqNo.dat', STATUS='replace'); write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='lastSobol.dat', STATUS='replace');  write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='legitSobol.dat', STATUS='replace');  write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='lastParam.dat', STATUS='replace');  write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='state.dat', STATUS='replace');  write(41,*) 0; close(41)
              open(UNIT=41, FILE='logFile.txt', STATUS='replace'); close(41)
              open(UNIT=41, FILE='internalConfig.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='sobol.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='sobolFnVal.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='missingSobol.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='x_starts.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='searchResults.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='searchStart.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='FinalResults.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='FinalStart.dat', STATUS='replace'); close(41)

            ELSEIF(OPTION==2) THEN
              !if updating Sobol points, then don't reset state, lastSobol, previously
              !calculated function values, or previously used parameters
              isUpdate= .TRUE.
              close(41)
              open(UNIT=41, FILE='missingSobol.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='seqNo.dat', STATUS='replace'); write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='lastParam.dat', STATUS='replace');  write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='x_starts.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='searchResults.dat', STATUS='replace');  close(41)
              open(UNIT=41, FILE='searchStart.dat', STATUS='replace');  close(41)
              open(UNIT=41, FILE='FinalResults.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='FinalStart.dat', STATUS='replace'); close(41)

              ! if we are updating, then store old values first
              call parseConfig(.FALSE., 'internalConfig.dat')
              old_p_qr_ndraw = p_qr_ndraw

            ELSEIF(OPTION==3) THEN
              !if updating number of local points, only reset  the following files
              isUpdate= .TRUE.
              close(41)
              open(UNIT=41, FILE='seqNo.dat', STATUS='replace'); write(41,*) 0 ; close(41)
              open(UNIT=41, FILE='x_starts.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='FinalResults.dat', STATUS='replace'); close(41)
              open(UNIT=41, FILE='FinalStart.dat', STATUS='replace'); close(41)

              ! if we are updating, then store old values first
              call parseConfig(.FALSE., 'internalConfig.dat')
              old_p_qr_ndraw = p_qr_ndraw
              old_p_maxpoints = p_maxpoints
            END IF

            call parseConfig(isUpdate, config)

            !print to internal config file so we have values for warm starts.
            call myopen(UNIT=fileDesc, FILE='internalConfig.dat', STATUS='replace', &
             IOSTAT=openStat, ACTION='write',SHARE='DENYRW')
            write(fileDesc,*) p_nx
            write(fileDesc,*) p_nmom
            write(fileDesc,*) p_maxeval
            write(fileDesc,*) p_qr_ndraw
            write(fileDesc,*) p_legitimate
            write(fileDesc,*) p_fvalmax
            write(fileDesc,*) p_maxpoints-1
            write(fileDesc,*) p_searchType
            write(fileDesc,*) p_lotteryPoints
            DO i=1,p_nx
                write(fileDesc,*) p_range(i,1),p_range(i,2)
            END DO
            DO i=1,p_nx
                write(fileDesc,*) p_init(i)
            END DO
            DO i=1,p_nx
                write(fileDesc,*) p_bound(i,1),p_bound(i,2)
            END DO
            call myclose(fileDesc)

            IF (OPTION==2) THEN  ! Update sobol points
                ! Then decide which sobol point to start with or
                ! to move on to local minimizations.
                INQUIRE(FILE='sobolFnVal.dat', EXIST=THERE_EXISTS)
                IF(THERE_EXISTS) THEN
                  allocate(Old_sobolFnVal(2*old_p_qr_ndraw,p_nx+2))
                  call myread2(old_sobolFnVal,'sobolFnVal.dat',numsolved)
                  numsolved=maxval(INT(old_sobolFnVal(1:numsolved,1)))
                  deallocate(Old_sobolFnVal)
                  legitSobol=getNumber('legitSobol.dat')
                  print*, "new # sobol= ",p_qr_ndraw
                  print*, "old # sobol= ",old_p_qr_ndraw
                  print*, "solved # sobol= ",numsolved
                  print*, "legitimate Sobol #= ",legitSobol
                  IF (numsolved >= p_qr_ndraw .or. legitSobol >= p_legitimate) THEN
                    write(errorString,*) "Update sobol: numsolved > p_qr_ndraw or legitSobol > p_legitimate ", &
                    NEW_LINE('A'), old_p_qr_ndraw, p_qr_ndraw, legitSobol , p_legitimate, NEW_LINE('A'), &
                    " Change p_qr_ndraw to ", max(numsolved+1,p_qr_ndraw)," or p_legitimate to ", max(legitSobol+1,p_legitimate)
                    call exitState(errorString);
                  ENDIF
                  CALL setState(numsolved, 'lastSobol.dat')
                ELSE
                  CALL setState(0, 'lastSobol.dat')
                ENDIF
            ENDIF

            IF (OPTION==3) THEN  ! Update number of local searches

                ! Decide  which sobol point to start the local minimization with.

                if(old_p_qr_ndraw .NE. p_qr_ndraw) THEN
                  print*, "new # sobol= ",p_qr_ndraw
                  print*, "old # sobol= ",old_p_qr_ndraw
                  write(errorString,*) "Update LocalMinimizations: old_p_qr_ndraw .NE. p_qr_ndraws", &
                  NEW_LINE('A'), "Run with option=2, update sobol points first "
                  call exitState(errorString);
                ENDIF

                INQUIRE(FILE='sobolFnVal.dat', EXIST=THERE_EXISTS)
                IF(THERE_EXISTS) THEN
                  allocate(Old_sobolFnVal(2*old_p_qr_ndraw,p_nx+2))
                  call myread2(old_sobolFnVal,'sobolFnVal.dat',numsolved)
                  numsolved=maxval(INT(old_sobolFnVal(1:numsolved,1)))
                  deallocate(Old_sobolFnVal)
                  print*, "solved # sobol= ",numsolved
                  print*, "new # local mimizations = ",p_maxpoints
                  IF (numsolved < p_maxpoints) THEN
                    write(errorString,*) "Update LocalMinimizations: numsolved < p_maxpoints", &
                    NEW_LINE('A'), "Run with option=2, update sobol points first!"
                    call exitState(errorString);
                  ENDIF
                ELSE
                  write(errorString,*) "Update LocalMinimizations:  sobolFnVal.dat cannot be found!"
                  call exitState(errorString);
                ENDIF

                INQUIRE(FILE='searchResults.dat', EXIST=THERE_EXISTS)
                IF(THERE_EXISTS) THEN
                  ALLOCATE(Old_SearchResults(2*old_p_maxpoints,p_nx+4))
                  call myread2(Old_SearchResults,'searchResults.dat',numsolved)
                  numsolved=maxval(INT(Old_SearchResults(1:numsolved,2)))
                  deallocate(Old_SearchResults)
                  print*, "new # local mimizations = ",p_maxpoints
                  print*, "old # local mimizations = ",old_p_maxpoints
                  print*, "solved # local mimizations= ",numsolved
                  IF (numsolved > p_maxpoints) THEN
                    write(errorString,*) "Update LocalMinimizations: numsolved > p_maxpoints", &
                    NEW_LINE('A'), old_p_maxpoints, p_maxpoints, &
                    " Change p_maxpoints to >", numsolved+1
                    call exitState(errorString);
                  ENDIF
                  CALL setState(numsolved, 'lastParam.dat')
                ELSE
                  CALL setState(0, 'lastParam.dat')
                ENDIF
            ENDIF

            if(p_qr_ndraw==0 .or. OPTION==3) THEN
              CALL setState(6, 'state.dat')
            ELSE
              CALL setState(1, 'state.dat')
            endif
        ELSE
            !Otherwise, read from internal config file
            call parseConfig(.FALSE., 'internalConfig.dat')
        END IF
        seqNo = getNextNumber('seqNo.dat')
    END SUBROUTINE initialize

    SUBROUTINE parseConfig(isUpdate, configFile)
        !parse the configuration file
        LOGICAL, INTENT(IN) :: isUpdate
        CHARACTER(LEN=*), INTENT(IN) :: configFile
        INTEGER(I4B) :: openStat, i
        INTEGER(I4B) :: nx, nmom, maxeval, ninterppt,legitimate
        REAL(DP) :: fvalmax
        INTEGER(I4B) :: qr_ndraw,maxpoints
        INTEGER(I4B) :: searchType, lotteryPoints, fileDesc
        CHARACTER(LEN=100) :: line
        LOGICAL :: parsedLine

        call myopen(UNIT=fileDesc, FILE=configFile, STATUS='old', IOSTAT=openStat, &
        ACTION='read',SHARE='DENYRW')
        IF (openStat /= 0) THEN
            ! file doesn't exist. This is an error
            print *,"<genericParams.parseConfig()> Error: unable to open file:",configFile
            stop 0
        END IF

        !Parse the basic parameters in the first line
        parsedLine = .FALSE.
        DO WHILE (parsedLine .eqv. .FALSE.)
            read(fileDesc,'(A100)',END=10) line

            IF (line(1:1)=='!') THEN
                cycle
            END IF

            parsedLine = .TRUE.

            read(line,*) nx
            read(fileDesc,*,END=10) nmom
            read(fileDesc,*,END=10) maxeval
            read(fileDesc,*,END=10) qr_ndraw
            read(fileDesc,*,END=10) legitimate
            read(fileDesc,*,END=10) fvalMax
            read(fileDesc,*,END=10) maxpoints
            read(fileDesc,*,END=10) searchType
            read(fileDesc,*,END=10) lotteryPoints
        END DO

        IF (maxeval < 0) THEN
            maxeval = 40*(nx+1)
        END IF

        ninterppt = 2*nx+1

        IF (qr_ndraw < 0) THEN
            qr_ndraw = 500
        END IF

        IF (legitimate < 0) THEN ! minimum # of legitimate draws is not binding.
            legitimate = qr_ndraw + 10000
        END IF

        IF (fvalMax < 0.0_DP) THEN ! All objective values are legitimate, no matter how big.
            fvalMax = 100000000.0_DP
        END IF

        IF (qr_ndraw == 0) THEN
            maxpoints = 0

        ELSEIF (maxpoints < 0 ) THEN
            maxpoints=min(200,qr_ndraw)

        ELSE IF(maxpoints > qr_ndraw) THEN
            print *, "<genericParams:parseConfig> Error in config file. maxpoints > # sobol points"
            stop 0
        END IF

        IF (lotteryPoints < 0) THEN
            lotteryPoints = maxpoints
        END IF

        IF (isUpdate) THEN
            !if updating, don't allow changes in nx or nmom
            IF (nx /= p_nx) THEN
                print *, "<genericParams:parseConfig> Cannot update the number of parameters.",nx,p_nx
                stop 0
            END IF

            IF (nmom /= p_nmom) THEN
                print *, "<genericParams:parseConfig> Cannot update the number of moments."
                stop 0
            END IF
        ELSE
            IF (allocated(p_range) .eqv. .FALSE.) THEN
                allocate(p_range(nx,2))
            END IF
            IF (allocated(p_init) .eqv. .FALSE.) THEN
                allocate(p_init(nx))
            END IF
            IF (allocated(p_bound) .eqv. .FALSE.) THEN
                allocate(p_bound(nx,2))
            END IF
        END IF

        p_nx=nx
        p_nmom=nmom
        p_maxeval=maxeval
        p_qr_ndraw=qr_ndraw
        p_legitimate = legitimate
        p_fvalmax = fvalMax
        p_maxpoints = maxpoints+1 !we add one to account for the initial guess
        p_ninterppt = ninterppt
        p_searchType = searchType
        p_lotteryPoints = lotteryPoints

        ! now parse the range for the sobol points
        parsedLine = .FALSE.
        DO WHILE ( (parsedLine .eqv. .FALSE.))
            read(fileDesc,'(A100)',END=11) line

            IF (line(1:1)=='!') THEN
                cycle
            END IF

            parsedLine = .TRUE.

            ! We are not allowed to change the parameter ranges except for warm start
            IF(.not. isUpdate) read(line,*) p_range(1,1), p_range(1,2)

            DO i=2,p_nx
                IF(.not. isUpdate) THEN
                read(fileDesc,*,END=11) p_range(i,1), p_range(i,2)
                ELSE
                  read(fileDesc,'(A100)',END=11) line
                ENDIF
            END DO
        END DO

        ! now parse the initial guess for the parameters
        ! initial guess can be updated during cold starts phase.
        parsedLine = .FALSE.
        DO WHILE ( (parsedLine .eqv. .FALSE.))
            read(fileDesc,'(A100)',END=12) line

            IF (line(1:1)=='!') THEN
                cycle
            END IF

            parsedLine = .TRUE.
            ! We will always change the initial parameter values
            read(line,*) p_init(1)

            DO i=2,p_nx
                read(fileDesc,*,END=12) p_init(i)
            END DO
        END DO

        ! now parse the bounds for the parameters
        parsedLine = .FALSE.
        DO WHILE ( (parsedLine .eqv. .FALSE.))
            read(fileDesc,'(A100)',END=13) line

            IF (line(1:1)=='!') THEN
                cycle
            END IF

            parsedLine = .TRUE.
            IF(.not. isUpdate)  read(line,*) p_bound(1,1), p_bound(1,2)

            DO i=2,p_nx
              IF(.not. isUpdate) THEN
                read(fileDesc,*,END=13) p_bound(i,1), p_bound(i,2)
              ELSE
                read(fileDesc,'(A100)',END=11) line
              ENDIF
            END DO
        END DO

        call myclose(fileDesc)

        IF(allocated(p_wspace)) deallocate(p_wspace)
        allocate(p_wspace((p_ninterppt+5)*(p_ninterppt+p_nx)+3*p_nx*(p_nx+5)/2))

        return

10      print *,"<genericParams:parseConfig> Error. ",configFile," Unable to read first line of parameters"
11      print *,"<genericParams:parseConfig> Error. ",configFile," Unable to read range of values for all parameters"
12      print *,"<genericParams:parseConfig> Error. ",configFile," Unable to read initial guesses for all parameters"
13      print *,"<genericParams:parseConfig> Error. ",configFile," Unable to read bounds for all parameters"
!100     format(i10)

    END SUBROUTINE parseConfig
END MODULE genericParams

max_stars_repo_name: serdarozkan/TikTak

max_issues_repo_name: serdarozkan/TikTak

max_forks_repo_name: serdarozkan/TikTak


!*==ssptrd.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b SSPTRD
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download SSPTRD + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE SSPTRD( UPLO, N, AP, D, E, TAU, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          UPLO
!       INTEGER            INFO, N
!       ..
!       .. Array Arguments ..
!       REAL               AP( * ), D( * ), E( * ), TAU( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> SSPTRD reduces a real symmetric matrix A stored in packed form to
!> symmetric tridiagonal form T by an orthogonal similarity
!> transformation: Q**T * A * Q = T.
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] UPLO
!> \verbatim
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!> \endverbatim
!>
!> \param[in,out] AP
!> \verbatim
!>          AP is REAL array, dimension (N*(N+1)/2)
!>          On entry, the upper or lower triangle of the symmetric matrix
!>          A, packed columnwise in a linear array.  The j-th column of A
!>          is stored in the array AP as follows:
!>          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
!>          if UPLO = 'L', AP(i + (j-1)*(2*n-j)/2) = A(i,j) for j<=i<=n.
!>          On exit, if UPLO = 'U', the diagonal and first superdiagonal
!>          of A are overwritten by the corresponding elements of the
!>          tridiagonal matrix T, and the elements above the first
!>          superdiagonal, with the array TAU, represent the orthogonal
!>          matrix Q as a product of elementary reflectors; if UPLO
!>          = 'L', the diagonal and first subdiagonal of A are over-
!>          written by the corresponding elements of the tridiagonal
!>          matrix T, and the elements below the first subdiagonal, with
!>          the array TAU, represent the orthogonal matrix Q as a product
!>          of elementary reflectors. See Further Details.
!> \endverbatim
!>
!> \param[out] D
!> \verbatim
!>          D is REAL array, dimension (N)
!>          The diagonal elements of the tridiagonal matrix T:
!>          D(i) = A(i,i).
!> \endverbatim
!>
!> \param[out] E
!> \verbatim
!>          E is REAL array, dimension (N-1)
!>          The off-diagonal elements of the tridiagonal matrix T:
!>          E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'.
!> \endverbatim
!>
!> \param[out] TAU
!> \verbatim
!>          TAU is REAL array, dimension (N-1)
!>          The scalar factors of the elementary reflectors (see Further
!>          Details).
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup realOTHERcomputational
!
!> \par Further Details:
!  =====================
!>
!> \verbatim
!>
!>  If UPLO = 'U', the matrix Q is represented as a product of elementary
!>  reflectors
!>
!>     Q = H(n-1) . . . H(2) H(1).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v**T
!>
!>  where tau is a real scalar, and v is a real vector with
!>  v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in AP,
!>  overwriting A(1:i-1,i+1), and tau is stored in TAU(i).
!>
!>  If UPLO = 'L', the matrix Q is represented as a product of elementary
!>  reflectors
!>
!>     Q = H(1) H(2) . . . H(n-1).
!>
!>  Each H(i) has the form
!>
!>     H(i) = I - tau * v * v**T
!>
!>  where tau is a real scalar, and v is a real vector with
!>  v(1:i) = 0 and v(i+1) = 1; v(i+2:n) is stored on exit in AP,
!>  overwriting A(i+2:n,i), and tau is stored in TAU(i).
!> \endverbatim
!>
!  =====================================================================
      SUBROUTINE SSPTRD(Uplo,N,Ap,D,E,Tau,Info)
      IMPLICIT NONE
!*--SSPTRD154
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Uplo
      INTEGER Info , N
!     ..
!     .. Array Arguments ..
      REAL Ap(*) , D(*) , E(*) , Tau(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      REAL ONE , ZERO , HALF
      PARAMETER (ONE=1.0,ZERO=0.0,HALF=1.0/2.0)
!     ..
!     .. Local Scalars ..
      LOGICAL upper
      INTEGER i , i1 , i1i1 , ii
      REAL alpha , taui
!     ..
!     .. External Subroutines ..
      EXTERNAL SAXPY , SLARFG , SSPMV , SSPR2 , XERBLA
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      REAL SDOT
      EXTERNAL LSAME , SDOT
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters
!
      Info = 0
      upper = LSAME(Uplo,'U')
      IF ( .NOT.upper .AND. .NOT.LSAME(Uplo,'L') ) THEN
         Info = -1
      ELSEIF ( N<0 ) THEN
         Info = -2
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('SSPTRD',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N<=0 ) RETURN
!
      IF ( upper ) THEN
!
!        Reduce the upper triangle of A.
!        I1 is the index in AP of A(1,I+1).
!
         i1 = N*(N-1)/2 + 1
         DO i = N - 1 , 1 , -1
!
!           Generate elementary reflector H(i) = I - tau * v * v**T
!           to annihilate A(1:i-1,i+1)
!
            CALL SLARFG(i,Ap(i1+i-1),Ap(i1),1,taui)
            E(i) = Ap(i1+i-1)
!
            IF ( taui/=ZERO ) THEN
!
!              Apply H(i) from both sides to A(1:i,1:i)
!
               Ap(i1+i-1) = ONE
!
!              Compute  y := tau * A * v  storing y in TAU(1:i)
!
               CALL SSPMV(Uplo,i,taui,Ap,Ap(i1),1,ZERO,Tau,1)
!
!              Compute  w := y - 1/2 * tau * (y**T *v) * v
!
               alpha = -HALF*taui*SDOT(i,Tau,1,Ap(i1),1)
               CALL SAXPY(i,alpha,Ap(i1),1,Tau,1)
!
!              Apply the transformation as a rank-2 update:
!                 A := A - v * w**T - w * v**T
!
               CALL SSPR2(Uplo,i,-ONE,Ap(i1),1,Tau,1,Ap)
!
               Ap(i1+i-1) = E(i)
            ENDIF
            D(i+1) = Ap(i1+i)
            Tau(i) = taui
            i1 = i1 - i
         ENDDO
         D(1) = Ap(1)
      ELSE
!
!        Reduce the lower triangle of A. II is the index in AP of
!        A(i,i) and I1I1 is the index of A(i+1,i+1).
!
         ii = 1
         DO i = 1 , N - 1
            i1i1 = ii + N - i + 1
!
!           Generate elementary reflector H(i) = I - tau * v * v**T
!           to annihilate A(i+2:n,i)
!
            CALL SLARFG(N-i,Ap(ii+1),Ap(ii+2),1,taui)
            E(i) = Ap(ii+1)
!
            IF ( taui/=ZERO ) THEN
!
!              Apply H(i) from both sides to A(i+1:n,i+1:n)
!
               Ap(ii+1) = ONE
!
!              Compute  y := tau * A * v  storing y in TAU(i:n-1)
!
               CALL SSPMV(Uplo,N-i,taui,Ap(i1i1),Ap(ii+1),1,ZERO,Tau(i),&
     &                    1)
!
!              Compute  w := y - 1/2 * tau * (y**T *v) * v
!
               alpha = -HALF*taui*SDOT(N-i,Tau(i),1,Ap(ii+1),1)
               CALL SAXPY(N-i,alpha,Ap(ii+1),1,Tau(i),1)
!
!              Apply the transformation as a rank-2 update:
!                 A := A - v * w**T - w * v**T
!
               CALL SSPR2(Uplo,N-i,-ONE,Ap(ii+1),1,Tau(i),1,Ap(i1i1))
!
               Ap(ii+1) = E(i)
            ENDIF
            D(i) = Ap(ii)
            Tau(i) = taui
            ii = i1i1
         ENDDO
         D(N) = Ap(ii)
      ENDIF
!
!
!     End of SSPTRD
!
      END SUBROUTINE SSPTRD


text:
Davis is wellknown for its fantastic environment for kids. The city is full of good places to take kids, plenty of which are fun for parents too.

 Within Davis

  Davis has a plethora of fantastic Parks
  Explorit Science Center
  Davis Farmers Market (and Picnic In The Park in the summer time), and the Davis Flea Market
  XCentral Park Play FountainX (currently turned off)
  Impossible Acres has farm animals, pumpkin patches before Halloween. During summer months it has a UPick patch for berries, cherries, apricots, etc.
  The The Flying Carousel of the Delta Breeze carousel in Central Park during Picnic in the Park and the Farmers Market.
  XThird and BX Unless your kids like looking at bicycles.
  The Davis Food CoOp is a great place to shop with kids.  They have minishopping carts and a coloring area.
  The MU Bowling Alley has bumpers, light balls and even a ramp suitable for even the youngest of kids.
  Davis has multiple public Swimming pools
  For a very simple Davis outing, the intersection of 2nd Street 2nd and H Street H streets features a PastaQ Fountain great fountain, Train Station trains, and Village Bakery pizza by the slice near lots of benches.  You cant lose.
  If you have a group of kids you can arrange a UC Davis Barn Tours UC Davis Barn Tour
  California Raptor Center has lots of cool birds, including a bald eagle. There is a oneroom educational museum packed with interesting activities like identifying bones pulled from owl pellets.

 Outside of Davis
(sorted by rough distance from Davis)

  Vic Fazio Yolo Wildlife Area
  Old Sacramento
   http://www.csrmf.org/ California State Railroad Museum  $35 yearround pass for a kid gets 1 adult in free with the kid. Conveniently located within walking distance of the Amtrak station for even more train fun. Appropriate for toddlers and up  children too young to enjoy the exhibits will probably have fun with the train table.
   http://www.thediscovery.org/ Gold Rush History Center  Very well set up for kids; Explorit Science Center Davis Explorit members get in free.
  The Sacramento Zoo  An excellent zoo well set up for kids (of course). http://www.funderlandpark.com/ Funderland and http://www.sacwiki.org/Fairytale_Town Fairytale Town are right next door.
  The Jelly Belly Factory in Fairfield gives free tours  Weekdays are best since thats when the cool machines are running.
  The http://www.sackids.org/ Sacramento Childrens Museum in Rancho Cordova is designed for ages birth to eightyearsold.
  The Nimbus Salmon and Steelhead Hatchery  Any time of year, but early November is the best

 Other Activities For Kids in Davis

  Kids love Geocaching
  Youth Sports

 Seasonal Events For Kids

 || April || Picnic Day has more things to do than you can fit into a week let alone a day. ||
 || May || Celebrate Davis hosts tons of stuff for kids. ||
 || 4th of July July 4th || Many events for kids, including a 4th of July/Kids Parade kids parade (2 blocks long) and, of course, July 4th Fireworks||
 || October 31 || Downtown turns into a wonderful shoptoshop trickortreat fest. ||
 || December || A candlewalk 2 blocks long to the Holiday Tree Lighting Ceremony ||

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program pg17
    implicit none
    integer, parameter :: mml=8
    double precision   :: mat_a(mml,mml)
    character(len=64)  :: dfmt
    write(dfmt, '(a,i0,a)') '(' , mml , 'f7.1)'

    mat_a        = 10
    mat_a(1,1)   = 0
    write(*,dfmt) transpose(mat_a)
end program pg17

max_stars_repo_name: rsoftone/fortran-training

max_issues_repo_name: rsoftone/fortran-training

max_forks_repo_name: rsoftone/fortran-training


text: Welcome to the Wiki Hola! Ni hao! Hi!