vdw_ts.F90

Go to the documentation of this file.

!! Copyright (C) 2015 X. Andrade, R. A. DiStasio Jr., B. Santra
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module vdw_ts_oct_m
  use atom_oct_m
  use debug_oct_m
  use global_oct_m
  use hirshfeld_oct_m
  use io_oct_m
  use io_function_oct_m
  use ions_oct_m
  use, intrinsic :: iso_fortran_env
  use lattice_vectors_oct_m
  use lalg_basic_oct_m
  use messages_oct_m
  use mesh_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use ps_oct_m
  use space_oct_m
  use species_oct_m
  use pseudopotential_oct_m
  use unit_oct_m
  use unit_system_oct_m
 
  implicit none
 
  private
 
  public ::                               &
    vdw_ts_t,                             &
    vdw_ts_init,                          &
    vdw_ts_end,                           &
    vdw_ts_write_c6ab,                    &
    vdw_ts_force_calculate,               &
    vdw_ts_calculate
 
  type vdw_ts_t
    private
    real(real64), allocatable :: c6free(:)
    real(real64), allocatable :: dpfree(:)
    real(real64), allocatable :: r0free(:)
    real(real64), allocatable :: c6abfree(:, :)
    real(real64), allocatable :: volfree(:)
    real(real64), allocatable :: c6ab(:, :)
    real(real64)       :: cutoff
    real(real64)       :: damping
    real(real64)       :: sr
 
    real(real64), allocatable :: derivative_coeff(:)
  end type vdw_ts_t
 
contains
 
  subroutine vdw_ts_init(this, namespace, ions)
    type(vdw_ts_t),      intent(out)   :: this
    type(namespace_t),   intent(in)    :: namespace
    type(ions_t),        intent(in)    :: ions
 
    integer :: ispecies, jspecies
    real(real64) :: num, den
 
    push_sub(vdw_ts_init)
 
    !%Variable VDW_TS_cutoff
    !%Type float
    !%Default 10.0
    !%Section Hamiltonian::XC
    !%Description
    !% Set the value of the cutoff (unit of length) for the VDW correction in periodic system
    !% in the Tkatchenko and Scheffler (vdw_ts) scheme only.
    !%End
    call parse_variable(namespace, 'VDW_TS_cutoff', 10.0_real64, this%cutoff, units_inp%length)
 
 
    !%Variable VDW_TS_damping
    !%Type float
    !%Default 20.0
    !%Section Hamiltonian::XC
    !%Description
    !% Set the value of the damping function (in unit of 1/length) steepness for the VDW correction in the
    !% Tkatchenko-Scheffler scheme. See Equation (12) of Phys. Rev. Lett. 102 073005 (2009).
    !%End
    call parse_variable(namespace, 'VDW_TS_damping', 20.0_real64, this%damping, units_inp%length**(-1))
 
    !%Variable VDW_TS_sr
    !%Type float
    !%Default 0.94
    !%Section Hamiltonian::XC
    !%Description
    !% Set the value of the sr parameter in the damping function of the VDW correction in the
    !% Tkatchenko-Scheffler scheme. See Equation (12) of Phys. Rev. Lett. 102 073005 (2009).
    !% This parameter depends on the xc functional used.
    !% The default value is 0.94, which holds for PBE. For PBE0, a value of 0.96 should be used.
    !%End
    call parse_variable(namespace, 'VDW_TS_sr', 0.94_real64, this%sr)
 
 
    safe_allocate(this%c6free(1:ions%nspecies))
    safe_allocate(this%dpfree(1:ions%nspecies))
    safe_allocate(this%r0free(1:ions%nspecies))
    safe_allocate(this%volfree(1:ions%nspecies))
    safe_allocate(this%c6abfree(1:ions%nspecies, 1:ions%nspecies))
    safe_allocate(this%c6ab(1:ions%natoms, 1:ions%natoms))
    safe_allocate(this%derivative_coeff(1:ions%natoms))
 
    do ispecies = 1, ions%nspecies
      call get_vdw_param(namespace, ions%species(ispecies)%s%get_label(), &
        this%c6free(ispecies), this%dpfree(ispecies), this%r0free(ispecies))
      select type(spec=>ions%species(ispecies)%s)
      type is(pseudopotential_t)
        this%volfree(ispecies) = ps_density_volume(spec%ps, namespace)
      class default
        assert(.false.)
      end select
    end do
 
    do ispecies = 1, ions%nspecies
      do jspecies = 1, ions%nspecies
        num = m_two*this%c6free(ispecies)*this%c6free(jspecies)
        den = (this%dpfree(jspecies)/this%dpfree(ispecies))*this%c6free(ispecies) &
          + (this%dpfree(ispecies)/this%dpfree(jspecies))*this%c6free(jspecies)
        this%c6abfree(ispecies, jspecies) = num/den
      end do
    end do
    pop_sub(vdw_ts_init)
  end subroutine vdw_ts_init
 
  !------------------------------------------
 
  subroutine vdw_ts_end(this)
    type(vdw_ts_t), intent(inout) :: this
 
    push_sub(vdw_ts_end)
 
    safe_deallocate_a(this%c6free)
    safe_deallocate_a(this%dpfree)
    safe_deallocate_a(this%r0free)
    safe_deallocate_a(this%volfree)
    safe_deallocate_a(this%c6abfree)
    safe_deallocate_a(this%c6ab)
    safe_deallocate_a(this%derivative_coeff)
 
    pop_sub(vdw_ts_end)
  end subroutine vdw_ts_end
 
  !------------------------------------------
  subroutine vdw_ts_calculate(this, namespace, space, latt, atom, natoms, pos, mesh, nspin, density, energy, potential, force)
    type(vdw_ts_t),           intent(inout) :: this
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    type(lattice_vectors_t),  intent(in)    :: latt
    type(atom_t),             intent(in)    :: atom(:)
    integer,                  intent(in)    :: natoms
    real(real64),             intent(in)    :: pos(1:space%dim,1:natoms)
    class(mesh_t),            intent(in)    :: mesh
    integer,                  intent(in)    :: nspin
    real(real64),             intent(in)    :: density(:, :)
    real(real64),             intent(out)   :: energy
    real(real64), contiguous, intent(out)   :: potential(:)
    real(real64), contiguous, intent(out)   :: force(:, :)
 
    interface
      subroutine f90_vdw_calculate(natoms, dd, sr, zatom, coordinates, vol_ratio, &
        energy, force, derivative_coeff)
        use, intrinsic :: iso_fortran_env
        implicit none
        integer, intent(in)  :: natoms
        real(real64), intent(in)  :: dd
        real(real64), intent(in)  :: sr
        integer,      intent(in)  :: zatom
        real(real64), intent(in)  :: coordinates
        real(real64), intent(in)  :: vol_ratio
        real(real64), intent(out) :: energy
        real(real64), intent(out) :: force
        real(real64), intent(out) :: derivative_coeff
      end subroutine f90_vdw_calculate
    end interface
 
    type(lattice_iterator_t) :: latt_iter
    integer :: iatom, jatom, ispecies, jspecies, jcopy, ip
    real(real64) :: rr, rr6, dffdr0, ee, ff, dee, dffdrab, dffdvra, deabdvra
    real(real64), allocatable :: coordinates(:,:), vol_ratio(:), dvadens(:), &
      dr0dvra(:), r0ab(:,:)
    type(hirshfeld_t) :: hirshfeld
    integer, allocatable :: zatom(:)
    real(real64) :: x_j(space%dim)
 
    push_sub(vdw_ts_calculate)
 
    safe_allocate(vol_ratio(1:natoms))
    safe_allocate(dvadens(1:mesh%np))
    safe_allocate(dr0dvra(1:natoms))
 
    energy=m_zero
    force(1:space%dim, 1:natoms) = m_zero
    this%derivative_coeff(1:natoms) = m_zero
    call hirshfeld_init(hirshfeld, mesh, namespace, space, latt, atom, natoms, pos, nspin)
 
    do iatom = 1, natoms
      call hirshfeld_volume_ratio(hirshfeld, iatom, density, vol_ratio(iatom))
    end do
 
    do iatom = 1, natoms
      ispecies = atom(iatom)%species%get_index()
      dr0dvra(iatom) = this%r0free(ispecies)/(m_three*(vol_ratio(iatom)**(m_two/m_three)))
      do jatom = 1, natoms
        jspecies = atom(jatom)%species%get_index()
        !the next operation is done again inside the .c part for the non periodic case
        this%c6ab(iatom,jatom) = vol_ratio(iatom)*vol_ratio(jatom)*this%c6abfree(ispecies,jspecies)
      end do
    end do
 
    if (space%is_periodic()) then ! periodic case
      safe_allocate(r0ab(1:natoms,1:natoms))
 
      !Precomputing some quantities
      do iatom = 1, natoms
        ispecies = atom(iatom)%species%get_index()
        do jatom = 1, natoms
          jspecies = atom(jatom)%species%get_index()
 
          r0ab(iatom,jatom) = (vol_ratio(iatom)**(m_one/m_three))*this%r0free(ispecies) &
            + (vol_ratio(jatom)**(m_one/m_three))*this%r0free(jspecies)
        end do
      end do
 
      latt_iter = lattice_iterator_t(latt, this%cutoff)
      do jatom = 1, natoms
        jspecies = atom(jatom)%species%get_index()
 
        do jcopy = 1, latt_iter%n_cells ! one of the periodic copy is the initial atom
          x_j = pos(:, jatom) + latt_iter%get(jcopy)
 
          do iatom = 1, natoms
            ispecies = atom(iatom)%species%get_index()
            rr =  norm2(x_j - pos(:, iatom))
            rr6 = rr**6
 
            if (rr < r_min_atom_dist) cycle !To avoid self interaction
 
            ee = exp(-this%damping * ((rr / (this%sr*r0ab(iatom, jatom))) - m_one))
            ff = m_one/(m_one + ee)
            dee = ee*ff**2
 
            !Calculate the derivative of the damping function with respect to the distance between atoms A and B.
            dffdrab = (this%damping/(this%sr*r0ab(iatom, jatom)))*dee
            !Calculate the derivative of the damping function with respect to the distance between the van der Waals radius.
            dffdr0 =  -this%damping*rr/(this%sr*r0ab(iatom, jatom)**2)*dee
 
            energy = energy - m_half*ff*this%c6ab(iatom, jatom)/rr6
 
            ! Derivative of the damping function with respecto to the volume ratio of atom A.
            dffdvra = dffdr0*dr0dvra(iatom)
 
            ! Calculation of the pair-wise partial energy derivative with respect to the volume ratio of atom A.
            deabdvra = (dffdvra*this%c6ab(iatom, jatom) + ff*vol_ratio(jatom)*this%c6abfree(ispecies, jspecies))/rr6
 
            this%derivative_coeff(iatom) = this%derivative_coeff(iatom) + deabdvra
 
          end do
        end do
      end do
 
      safe_deallocate_a(r0ab)
    else ! Non periodic case
      safe_allocate(coordinates(1:space%dim, 1:natoms))
      safe_allocate(zatom(1:natoms))
 
      do iatom = 1, natoms
        coordinates(:, iatom) = pos(:, iatom)
        zatom(iatom) = int(atom(iatom)%species%get_z())
 
      end do
 
      call f90_vdw_calculate(natoms,  this%damping, this%sr, zatom(1), coordinates(1, 1), &
        vol_ratio(1), energy, force(1, 1), this%derivative_coeff(1))
 
 
      safe_deallocate_a(coordinates)
      safe_deallocate_a(zatom)
    end if
 
    ! Calculate the potential
    potential = m_zero
    do iatom = 1, natoms
      call hirshfeld_density_derivative(hirshfeld, iatom, dvadens)
      call lalg_axpy(mesh%np, -this%derivative_coeff(iatom), dvadens, potential)
    end do
 
    if (debug%info) then
      call dio_function_output(1_8, "./", "vvdw", namespace, space, mesh, potential, unit_one, ip)
    end if
 
    call hirshfeld_end(hirshfeld)
 
    safe_deallocate_a(vol_ratio)
    safe_deallocate_a(dvadens)
    safe_deallocate_a(dr0dvra)
 
    pop_sub(vdw_ts_calculate)
  end subroutine vdw_ts_calculate
 
 
 
  !------------------------------------------
  subroutine vdw_ts_force_calculate(this, force_vdw, ions, mesh, nspin, density)
    type(vdw_ts_t),      intent(in)    :: this
    type(ions_t),        intent(in)    :: ions
    real(real64),        intent(inout) :: force_vdw(1:ions%space%dim, 1:ions%natoms)
    class(mesh_t),       intent(in)    :: mesh
    integer,             intent(in)    :: nspin
    real(real64),        intent(in)    :: density(:, :)
 
    type(hirshfeld_t) :: hirshfeld
    type(lattice_iterator_t) :: latt_iter
 
    integer :: iatom, jatom, ispecies, jspecies, jcopy
    real(real64) :: rr, rr6,  dffdr0, ee, ff, dee, dffdvra, deabdvra, deabdrab, x_j(ions%space%dim)
    real(real64), allocatable ::  vol_ratio(:), dvadrr(:), dr0dvra(:), r0ab(:,:), derivative_coeff(:), c6ab(:,:)
 
    push_sub(vdw_ts_force_calculate)
 
    call profiling_in("FORCE_VDW_TS")
 
    safe_allocate(vol_ratio(1:ions%natoms))
    safe_allocate(dvadrr(1:3))
    safe_allocate(dr0dvra(1:ions%natoms))
    safe_allocate(r0ab(1:ions%natoms,1:ions%natoms))
    safe_allocate(derivative_coeff(1:ions%natoms))
    safe_allocate(c6ab(1:ions%natoms,1:ions%natoms))
 
 
    force_vdw(1:ions%space%dim, 1:ions%natoms) = m_zero
    derivative_coeff(1:ions%natoms) = m_zero
    c6ab(1:ions%natoms,1:ions%natoms) = m_zero
    r0ab(1:ions%natoms,1:ions%natoms) = m_zero
    dr0dvra(1:ions%natoms) = m_zero
    dvadrr(1:ions%space%dim) = m_zero
    vol_ratio(1:ions%natoms) = m_zero
 
 
    call hirshfeld_init(hirshfeld, mesh, ions%namespace, ions%space, ions%latt, ions%atom, ions%natoms, ions%pos, nspin)
 
 
    do iatom = 1, ions%natoms
      call hirshfeld_volume_ratio(hirshfeld, iatom, density, vol_ratio(iatom))
    end do
 
    do iatom = 1, ions%natoms
      ispecies = ions%atom(iatom)%species%get_index()
      dr0dvra(iatom) = this%r0free(ispecies)/(m_three*(vol_ratio(iatom)**(m_two/m_three)))
      do jatom = 1, ions%natoms
        jspecies = ions%atom(jatom)%species%get_index()
        c6ab(iatom, jatom) = vol_ratio(iatom)*vol_ratio(jatom)*this%c6abfree(ispecies, jspecies)
      end do
    end do
 
    !Precomputing some quantities
    do iatom = 1, ions%natoms
      ispecies = ions%atom(iatom)%species%get_index()
      do jatom = iatom, ions%natoms
        jspecies = ions%atom(jatom)%species%get_index()
 
        r0ab(iatom, jatom) = (vol_ratio(iatom)**(m_one/m_three))*this%r0free(ispecies) &
          + (vol_ratio(jatom)**(m_one/m_three))*this%r0free(jspecies)
        if (iatom /= jatom) r0ab(jatom, iatom) = r0ab(iatom, jatom)
      end do
    end do
 
    latt_iter = lattice_iterator_t(ions%latt, this%cutoff)
    do jatom = 1, ions%natoms
      jspecies = ions%atom(jatom)%species%get_index()
 
      do jcopy = 1, latt_iter%n_cells ! one of the periodic copy is the initial atom
        x_j = ions%pos(:, jatom) + latt_iter%get(jcopy)
        do iatom = 1, ions%natoms
          ispecies = ions%atom(iatom)%species%get_index()
          rr =  norm2(x_j - ions%pos(:, iatom))
          rr6 = rr**6
 
          if (rr < r_min_atom_dist) cycle !To avoid self interaction
 
          ee = exp(-this%damping*(rr/(this%sr*r0ab(iatom, jatom)) - m_one))
          ff = m_one/(m_one + ee)
          dee = ee*ff**2
          !Calculate the derivative of the damping function with respect to the van der Waals radius.
          dffdr0 =  -this%damping*rr/( this%sr*r0ab(iatom, jatom)**2)*dee
          ! Calculation of the pair-wise partial energy derivative with respect to the distance between atoms A and B.
          deabdrab = c6ab(iatom,jatom)*(this%damping/(this%sr*r0ab(iatom, jatom))*dee - 6.0_real64*ff/rr)/rr6
          ! Derivative of the damping function with respecto to the volume ratio of atom A.
          dffdvra = dffdr0*dr0dvra(iatom)
          ! Calculation of the pair-wise partial energy derivative with respect to the volume ratio of atom A.
          deabdvra = (dffdvra*c6ab(iatom, jatom) + ff*vol_ratio(jatom)*this%c6abfree(ispecies, jspecies))/rr6
          !Summing for using later
          derivative_coeff(iatom) = derivative_coeff(iatom) + deabdvra
 
          ! Calculation of the pair-wise partial energy derivative with respect to the distance between atoms A and B.
          force_vdw(:, iatom) = force_vdw(:, iatom) + m_half*deabdrab*(ions%pos(:, iatom) - x_j)/rr
        end do
      end do
    end do
 
    do iatom = 1, ions%natoms
      do jatom = 1, ions%natoms
        ! dvadrr_ij = \frac{\delta V_i}{\delta \vec{x_j}}
        call hirshfeld_position_derivative(hirshfeld, ions%space, iatom, jatom, density, dvadrr)
        ! ions%atom(jatom)%f_vdw = sum_i coeff_i * dvadrr_ij
        force_vdw(:, jatom)= force_vdw(:, jatom) + derivative_coeff(iatom)*dvadrr
      end do
    end do
 
    call hirshfeld_end(hirshfeld)
 
    safe_deallocate_a(vol_ratio)
    safe_deallocate_a(dvadrr)
    safe_deallocate_a(dr0dvra)
    safe_deallocate_a(r0ab)
    safe_deallocate_a(derivative_coeff)
    safe_deallocate_a(c6ab)
 
    call profiling_out("FORCE_VDW_TS")
 
    pop_sub(vdw_ts_force_calculate)
  end subroutine vdw_ts_force_calculate
 
  !------------------------------------------
  subroutine vdw_ts_write_c6ab(this, ions, dir, fname, namespace)
    type(vdw_ts_t),      intent(in) :: this
    type(ions_t),        intent(in) :: ions
    character(len=*),    intent(in) :: dir, fname
    type(namespace_t),   intent(in) :: namespace
 
    integer :: iunit, iatom, jatom
 
    push_sub(vdw_ts_write_c6ab)
 
    if (ions%grp%is_root()) then
      call io_mkdir(dir, namespace)
      iunit = io_open(trim(dir) // "/" // trim(fname), namespace, action='write')
      write(iunit, '(a)') ' # Atom1 Atom2 C6_{12}^{eff}'
 
 
      do iatom = 1, ions%natoms
        do jatom = 1, ions%natoms
          write(iunit, '(3x, i5, i5, e15.6)') iatom, jatom, this%c6ab(iatom, jatom)
        end do
      end do
      call io_close(iunit)
    end if
    pop_sub(vdw_ts_write_c6ab)
 
  end subroutine vdw_ts_write_c6ab
 
 
  !------------------------------------------
 
  subroutine get_vdw_param(namespace, atom, c6, alpha, r0)
    type(namespace_t), intent(in)  :: namespace
    character(len=*),  intent(in)  :: atom
    real(real64),      intent(out) :: c6
    real(real64),      intent(out) :: alpha
    real(real64),      intent(out) :: r0
 
    push_sub(get_vdw_param)
 
    select case (trim(atom))
 
    case ('H')
      alpha = 4.500000_real64
      c6 = 6.500000_real64
      r0 = 3.100000_real64
 
    case ('He')
      alpha = 1.380000_real64
      c6 = 1.460000_real64
      r0 = 2.650000_real64
 
    case ('Li')
      alpha = 164.200000_real64
      c6 = 1387.000000_real64
      r0 = 4.160000_real64
 
    case ('Be')
      alpha = 38.000000_real64
      c6 = 214.000000_real64
      r0 = 4.170000_real64
 
    case ('B')
      alpha = 21.000000_real64
      c6 = 99.500000_real64
      r0 = 3.890000_real64
 
    case ('C')
      alpha = 12.000000_real64
      c6 = 46.600000_real64
      r0 = 3.590000_real64
 
    case ('N')
      alpha = 7.400000_real64
      c6 = 24.200000_real64
      r0 = 3.340000_real64
 
    case ('O')
      alpha = 5.400000_real64
      c6 = 15.600000_real64
      r0 = 3.190000_real64
 
    case ('F')
      alpha = 3.800000_real64
      c6 = 9.520000_real64
      r0 = 3.040000_real64
 
    case ('Ne')
      alpha = 2.670000_real64
      c6 = 6.380000_real64
      r0 = 2.910000_real64
 
    case ('Na')
      alpha = 162.700000_real64
      c6 = 1556.000000_real64
      r0 = 3.730000_real64
 
    case ('Mg')
      alpha = 71.000000_real64
      c6 = 627.000000_real64
      r0 = 4.270000_real64
 
    case ('Al')
      alpha = 60.000000_real64
      c6 = 528.000000_real64
      r0 = 4.330000_real64
 
    case ('Si')
      alpha = 37.000000_real64
      c6 = 305.000000_real64
      r0 = 4.200000_real64
 
    case ('P')
      alpha = 25.000000_real64
      c6 = 185.000000_real64
      r0 = 4.010000_real64
 
    case ('S')
      alpha = 19.600000_real64
      c6 = 134.000000_real64
      r0 = 3.860000_real64
 
    case ('Cl')
      alpha = 15.000000_real64
      c6 = 94.600000_real64
      r0 = 3.710000_real64
 
    case ('Ar')
      alpha = 11.100000_real64
      c6 = 64.300000_real64
      r0 = 3.550000_real64
 
    case ('K')
      alpha = 292.900000_real64
      c6 = 3897.000000_real64
      r0 = 3.710000_real64
 
    case ('Ca')
      alpha = 160.000000_real64
      c6 = 2221.000000_real64
      r0 = 4.650000_real64
 
    case ('Sc')
      alpha = 120.000000_real64
      c6 = 1383.000000_real64
      r0 = 4.590000_real64
 
    case ('Ti')
      alpha = 98.000000_real64
      c6 = 1044.000000_real64
      r0 = 4.510000_real64
 
    case ('V')
      alpha = 84.000000_real64
      c6 = 832.000000_real64
      r0 = 4.440000_real64
 
    case ('Cr')
      alpha = 78.000000_real64
      c6 = 602.000000_real64
      r0 = 3.990000_real64
 
    case ('Mn')
      alpha = 63.000000_real64
      c6 = 552.000000_real64
      r0 = 3.970000_real64
 
    case ('Fe')
      alpha = 56.000000_real64
      c6 = 482.000000_real64
      r0 = 4.230000_real64
 
    case ('Co')
      alpha = 50.000000_real64
      c6 = 408.000000_real64
      r0 = 4.180000_real64
 
    case ('Ni')
      alpha = 48.000000_real64
      c6 = 373.000000_real64
      r0 = 3.820000_real64
 
    case ('Cu')
      alpha = 42.000000_real64
      c6 = 253.000000_real64
      r0 = 3.760000_real64
 
    case ('Zn')
      alpha = 40.000000_real64
      c6 = 284.000000_real64
      r0 = 4.020000_real64
 
    case ('Ga')
      alpha = 60.000000_real64
      c6 = 498.000000_real64
      r0 = 4.190000_real64
 
    case ('Ge')
      alpha = 41.000000_real64
      c6 = 354.000000_real64
      r0 = 4.200000_real64
 
    case ('As')
      alpha = 29.000000_real64
      c6 = 246.000000_real64
      r0 = 4.110000_real64
 
    case ('Se')
      alpha = 25.000000_real64
      c6 = 210.000000_real64
      r0 = 4.040000_real64
 
    case ('Br')
      alpha = 20.000000_real64
      c6 = 162.000000_real64
      r0 = 3.930000_real64
 
    case ('Kr')
      alpha = 16.800000_real64
      c6 = 129.600000_real64
      r0 = 3.820000_real64
 
    case ('Rb')
      alpha = 319.200000_real64
      c6 = 4691.000000_real64
      r0 = 3.720000_real64
 
    case ('Sr')
      alpha = 199.000000_real64
      c6 = 3170.000000_real64
      r0 = 4.540000_real64
 
    case ('Y')
      alpha = 126.7370_real64
      c6 = 1968.580_real64
      r0 = 4.81510_real64
 
    case ('Zr')
      alpha = 119.97_real64
      c6 = 1677.91_real64
      r0 = 4.53_real64
 
    case ('Nb')
      alpha = 101.603_real64
      c6 = 1263.61_real64
      r0 = 4.2365_real64
 
    case ('Mo')
      alpha = 88.4225785_real64
      c6 = 1028.73_real64
      r0 = 4.099_real64
 
    case ('Tc')
      alpha = 80.083_real64
      c6 = 1390.87_real64
      r0 = 4.076_real64
 
    case ('Ru')
      alpha = 65.8950_real64
      c6 = 609.754_real64
      r0 = 3.99530_real64
 
 
    case ('Rh')
      alpha = 56.1_real64
      c6 = 469.0_real64
      r0 = 3.95_real64
 
    case ('Pd')
      alpha = 23.680000_real64
      c6 = 157.500000_real64
      r0 = 3.66000_real64
 
    case ('Ag')
      alpha = 50.600000_real64
      c6 = 339.000000_real64
      r0 = 3.820000_real64
 
    case ('Cd')
      alpha = 39.7_real64
      c6 = 452.0_real64
      r0 = 3.99_real64
    case ('Sn')
      alpha = 55.9500_real64
      c6 = 587.41700_real64
      r0 = 4.303000_real64
      !
    case ('Sb')
      alpha = 43.671970_real64
      c6 = 459.322_real64
      r0 = 4.2760_real64
 
 
    case ('Te')
      alpha = 37.65_real64
      c6 = 396.0_real64
      r0 = 4.22_real64
 
    case ('I')
      alpha = 35.000000_real64
      c6 = 385.000000_real64
      r0 = 4.170000_real64
 
    case ('Xe')
      alpha = 27.300000_real64
      c6 = 285.900000_real64
      r0 = 4.080000_real64
 
    case ('Cs')
      alpha = 427.12_real64
      c6 = 6582.08_real64
      r0 = 3.78_real64
 
    case ('Ba')
      alpha = 275.0_real64
      c6 = 5727.0_real64
      r0 = 4.77_real64
 
    case ('Hf')
      alpha = 99.52_real64
      c6 = 1274.8_real64
      r0 = 4.21_real64
 
    case ('Ta')
      alpha = 82.53_real64
      c6 = 1019.92_real64
      r0 = 4.15_real64
 
    case ('W')
      alpha = 71.041_real64
      c6 = 847.93_real64
      r0 = 4.08_real64
 
    case ('Re')
      alpha = 63.04_real64
      c6 = 710.2_real64
      r0 = 4.02_real64
 
    case ('Os')
      alpha = 55.055_real64
      c6 = 596.67_real64
      r0 = 3.84_real64
 
    case ('Ir')
      alpha = 42.51_real64
      c6 = 359.1_real64
      r0 = 4.00_real64
 
    case ('Pt')
      alpha = 39.68_real64
      c6 = 347.1_real64
      r0 = 3.92_real64
 
    case ('Au')
      alpha = 36.5_real64
      c6 = 298.0_real64
      r0 = 3.86_real64
 
    case ('Hg')
      alpha = 33.9_real64
      c6 = 392.0_real64
      r0 = 3.98_real64
 
    case ('Tl')
      alpha = 69.92_real64
      c6 = 717.44_real64
      r0 = 3.91_real64
 
    case ('Pb')
      alpha = 61.8_real64
      c6 = 697.0_real64
      r0 = 4.31_real64
 
    case ('Bi')
      alpha = 49.02_real64
      c6 = 571.0_real64
      r0 = 4.32_real64
 
    case ('Po')
      alpha = 45.013_real64
      c6 = 530.92_real64
      r0 = 4.097_real64
 
    case ('At')
      alpha = 38.93_real64
      c6 = 457.53_real64
      r0 = 4.07_real64
 
    case ('Rn')
      alpha = 33.54_real64
      c6 = 390.63_real64
      r0 = 4.23_real64
 
 
    case default
 
      call messages_write('vdw ts: reference free atom parameters not available for species '//trim(atom))
      call messages_fatal(namespace=namespace)
 
    end select
 
    pop_sub(get_vdw_param)
  end subroutine get_vdw_param
 
end module vdw_ts_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: