prg_xlbo_mod.F90

Go to the documentation of this file.

 
module prg_xlbo_mod
 
  use prg_openfiles_mod
  use bml
  use prg_kernelparser_mod
 
  implicit none
 
  private
 
  integer, parameter :: dp = kind(1.0d0)
 
  real(dp), parameter :: c0 = -6.0_dp
  real(dp), parameter :: c1 = 14.0_dp
  real(dp), parameter :: c2 = -8.0_dp
  real(dp), parameter :: c3 = -3.0_dp
  real(dp), parameter :: c4 = 4.0_dp
  real(dp), parameter :: c5 = -1.0_dp;
 
  real(dp), parameter :: kappa = 1.82_dp;
  real(dp), parameter :: alpha = 0.018_dp;
  real(dp), parameter :: cc = 0.9_dp;        ! Scaled prg_delta kernel
 
  type, public :: xlbo_type
 
    character(20) :: jobname
 
    integer :: verbose
 
    integer :: maxscfiter
 
    integer :: maxscfinititer
 
    real(dp) :: threshold
 
    integer :: minit
 
    real(dp) :: cc
 
  end type xlbo_type
 
  public :: prg_parse_xlbo, prg_xlbo_nint, prg_xlbo_nint_kernel, prg_xlbo_fcoulupdate
  public :: prg_xlbo_nint_kerneltimesres
 
contains
 
  subroutine prg_parse_xlbo(xlbo,filename)
 
    implicit none
    type(xlbo_type), intent(inout) :: xlbo
    integer, parameter :: nkey_char = 1, nkey_int = 4, nkey_re = 2, nkey_log = 1
    character(len=*) :: filename
 
    !Library of keywords with the respective defaults.
    character(len=50), parameter :: keyvector_char(nkey_char) = [character(len=100) :: &
         'JobName=']
    character(len=100) :: valvector_char(nkey_char) = [character(len=100) :: &
         'MyJob']
 
    character(len=50), parameter :: keyvector_int(nkey_int) = [character(len=50) :: &
         'Verbose=','Mprg_init=','MaxSCFIter=', 'MaxSCFInitIter=']
    integer :: valvector_int(nkey_int) = (/ &
         0, 6, 0, 4 /)
 
    character(len=50), parameter :: keyvector_re(nkey_re) = [character(len=50) :: &
         'NumThresh=', 'ScaledDeltaKernel=' ]
    real(dp) :: valvector_re(nkey_re) = (/&
         0.0, 0.99 /)
 
    character(len=50), parameter :: keyvector_log(nkey_log) = [character(len=100) :: &
         'Log1=']
    logical :: valvector_log(nkey_log) = (/&
         .false. /)
 
    !Start and stop characters
    character(len=50), parameter :: startstop(2) = [character(len=50) :: &
         'XLBO{', '}']
 
    call prg_parsing_kernel(keyvector_char,valvector_char&
         ,keyvector_int,valvector_int,keyvector_re,valvector_re,&
         keyvector_log,valvector_log,trim(filename),startstop)
 
    !Characters
    xlbo%JobName = valvector_char(1)
 
    !Reals
    xlbo%threshold = valvector_re(1)
    xlbo%cc = valvector_re(2)
 
    !Logicals
 
    !Integers
    xlbo%verbose = valvector_int(1)
    xlbo%minit = valvector_int(2)
    xlbo%maxscfiter = valvector_int(3)
    xlbo%maxscfinititer = valvector_int(4)
 
  end subroutine prg_parse_xlbo
 
 
  subroutine prg_xlbo_nint(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,xl)
    implicit none
    real(dp), allocatable, intent(inout) :: n(:), n_0(:), n_1(:), n_2(:), n_3(:), n_4(:), n_5(:)
    real(dp), allocatable, intent(in) :: charges(:)
    type(xlbo_type), intent(in) :: xl
 
    integer, intent(in) :: mdstep
    integer :: nats
 
    nats = size(charges,dim=1)
 
    if(.not.allocated(n))then
      allocate(n(nats))
      allocate(n_0(nats))
      allocate(n_1(nats))
      allocate(n_2(nats))
      allocate(n_3(nats))
      allocate(n_4(nats))
      allocate(n_5(nats))
    endif
 
    if(mdstep.le.1)then
      n = charges;
      n_0 = charges;
      n_1 = charges;
      n_2 = charges;
      n_3 = charges;
      n_4 = charges;
      n_5 = charges;
    endif
 
    n = 2.0_dp*n_0 - n_1 + xl%cc*kappa*(charges-n) &
         + alpha*(c0*n_0+c1*n_1+c2*n_2+c3*n_3+c4*n_4+c5*n_5);
    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n;
 
  end subroutine prg_xlbo_nint
 
 
  subroutine prg_xlbo_nint_kernel(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernel,xl)
    implicit none
    real(dp), allocatable, intent(inout) :: n(:), n_0(:), n_1(:), n_2(:), n_3(:), n_4(:), n_5(:)
    real(dp), allocatable, intent(in) :: charges(:)
    real(dp), allocatable, intent(in) :: kernel(:,:)
    type(xlbo_type), intent(in) :: xl
 
    integer, intent(in) :: mdstep
    integer :: nats
 
    nats = size(charges,dim=1)
 
    if(.not.allocated(n))then
      allocate(n(nats))
      allocate(n_0(nats))
      allocate(n_1(nats))
      allocate(n_2(nats))
      allocate(n_3(nats))
      allocate(n_4(nats))
      allocate(n_5(nats))
    endif
 
    if(mdstep.le.1)then
      n = charges;
      n_0 = charges;
      n_1 = charges;
      n_2 = charges;
      n_3 = charges;
      n_4 = charges;
      n_5 = charges;
    endif
 
    ! From developper's code
    !   dn2dt2 = -MATMUL(KK0,(q-n))
    !   n = 2*n_0 - n_1 + kappa*dn2dt2 +
    !   alpha*(C0*n_0+C1*n_1+C2*n_2+C3*n_3+C4*n_4+C5*n_5+C6*n_6)
    !   n_6 = n_5; n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n
 
    !call bml_print_matrix("ker",kernel,1,10,1,10)
    !write(*,*)matmul(kernel,(charges-n))
    !n = 2.0_dp*n_0 - n_1 + xl%cc*kappa*(charges-n) &
    n = 2.0_dp*n_0 - n_1 - 1.0_dp*kappa*matmul(kernel,(charges-n)) &
         + alpha*(c0*n_0+c1*n_1+c2*n_2+c3*n_3+c4*n_4+c5*n_5);
    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n;
 
  end subroutine prg_xlbo_nint_kernel
 
 
  subroutine prg_xlbo_nint_kerneltimesres(charges,n,n_0,n_1,n_2,n_3,n_4,n_5,mdstep,kernelTimesRes,xl)
    implicit none
    real(dp), allocatable, intent(inout) :: n(:), n_0(:), n_1(:), n_2(:), n_3(:), n_4(:), n_5(:)
    real(dp), allocatable, intent(in) :: charges(:)
    real(dp), allocatable, intent(in) :: kerneltimesres(:)
    type(xlbo_type), intent(in) :: xl
 
    integer, intent(in) :: mdstep
    integer :: nats
 
    nats = size(charges,dim=1)
 
    if(.not.allocated(n))then
      allocate(n(nats))
      allocate(n_0(nats))
      allocate(n_1(nats))
      allocate(n_2(nats))
      allocate(n_3(nats))
      allocate(n_4(nats))
      allocate(n_5(nats))
    endif
 
    if(mdstep.le.1)then
      n = charges;
      n_0 = charges;
      n_1 = charges;
      n_2 = charges;
      n_3 = charges;
      n_4 = charges;
      n_5 = charges;
    endif
 
    n = 2.0_dp*n_0 - n_1 - 1.0_dp*kappa*kerneltimesres &
         & + alpha*(c0*n_0+c1*n_1+c2*n_2+c3*n_3+c4*n_4+c5*n_5);
    n_5 = n_4; n_4 = n_3; n_3 = n_2; n_2 = n_1; n_1 = n_0; n_0 = n;
 
  end subroutine prg_xlbo_nint_kerneltimesres
 
 
 
  subroutine prg_xlbo_fcoulupdate(fcoul,charges,n)
    implicit none
    real(dp), intent(inout) :: fcoul(:,:),charges(:)
    real(dp), intent(inout) :: n(:)
 
    fcoul(1,:) = (2.0_dp*charges(:)-n(:))*fcoul(1,:)/n(:);
    fcoul(2,:) = (2.0_dp*charges(:)-n(:))*fcoul(2,:)/n(:);
    fcoul(3,:) = (2.0_dp*charges(:)-n(:))*fcoul(3,:)/n(:);
 
  end subroutine prg_xlbo_fcoulupdate
 
end module prg_xlbo_mod