modelmb_particles.F90

Go to the documentation of this file.

!! Copyright (C) 2009 N. Helbig and M. Verstraete
!!
!! 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 modelmb_particles_oct_m
  use, intrinsic :: iso_fortran_env
  use debug_oct_m
  use global_oct_m
  use messages_oct_m
  use namespace_oct_m
  use parser_oct_m
  use profiling_oct_m
  use space_oct_m
 
  implicit none
 
  private
 
  public ::                      &
    modelmb_particles_init,      &
    modelmb_particles_end,       &
    modelmb_particles_copy,      &
    modelmb_copy_masses,         &
    modelmb_particle_t
 
  type modelmb_particle_t
    private
    integer, public :: ndim
 
    integer, public :: ntype_of_particle
    integer :: max_particles_per_type
 
    integer, public :: nparticle = 0
 
    integer :: ndensities_to_calculate
 
    character(80), allocatable :: labels_particles(:)
 
    integer, allocatable, public :: particletype(:)
    integer, allocatable, public :: nparticles_per_type(:)
    integer, allocatable, public :: particles_of_type(:,:)
    integer, allocatable, public :: bosonfermion(:)
 
    integer, allocatable :: exchange_symmetry(:,:,:)
 
    real(real64), allocatable :: mass_particle(:)
 
    real(real64), allocatable, public :: charge_particle(:)
 
    character(80), allocatable :: labels_densities(:)
 
    integer, allocatable :: particle_kept_densities(:)
 
    integer,      allocatable, public :: nspindown(:,:)
    !                                                      for each type and state
    integer,      allocatable, public :: iyoung(:,:)
    real(real64), allocatable, public :: proj(:)
 
  contains
    procedure :: copy_masses => modelmb_copy_masses
  end type modelmb_particle_t
 
contains
 
  subroutine modelmb_particles_init(this, namespace, space, nst)
    type(modelmb_particle_t), intent(inout) :: this
    type(namespace_t),        intent(in)    :: namespace
    class(space_t),           intent(in)    :: space
    integer,                  intent(in)    :: nst
 
    integer :: ipart, ncols, nline, itmp, jtmp, ntype
    type(block_t) :: blk
 
    push_sub(modelmb_particles_init)
 
    ! read in scalar dimensions
 
    !%Variable NParticleModelmb
    !%Type integer
    !%Section States::ModelMB
    !%Default  0
    !%Description
    !% Number of particles in modelmb space.
    !% Full Ndim = <tt>NDimModelmb</tt>*<tt>NParticleModelmb</tt>
    !%End
    call parse_variable(namespace, 'NParticleModelmb', 0, this%nparticle)
 
    if (this%nparticle == 0) then
      pop_sub(modelmb_particles_init)
      return
    end if
 
    call messages_print_var_value("NParticleModelmb", this%nparticle, namespace=namespace)
 
    !%Variable NDimModelmb
    !%Type integer
    !%Section States::ModelMB
    !%Default 1
    !%Description
    !% Number of dimensions for modelmb space.
    !% Full Ndim = <tt>NDimModelmb</tt>*<tt>NParticleModelmb</tt>
    !%
    !%End
    call parse_variable(namespace, 'NDimModelmb', 1, this%ndim)
    call messages_print_var_value("NDimModelmb", this%ndim, namespace=namespace)
 
    !%Variable NTypeParticleModelmb
    !%Type integer
    !%Section States::ModelMB
    !%Default 1
    !%Description
    !% Number of different types of particles in modelmb space.
    !%End
    call parse_variable(namespace, 'NTypeParticleModelmb', 1, this%ntype_of_particle)
    call messages_print_var_value("NTypeParticleModelmb", this%ntype_of_particle, namespace=namespace)
    if (this%ntype_of_particle > this%nparticle) then
      write (message(1), '(2a,2I6)') ' Number of types of modelmb particles should be <= Number of modelmb particles ', &
        this%ntype_of_particle, this%nparticle
      call messages_fatal(1, namespace=namespace)
    end if
 
    if (this%ndim*this%nparticle /= space%dim) then
      message(1) = ' Number of modelmb particles * dimension of modelmb space must be = Ndim'
      call messages_fatal(1, namespace=namespace)
    end if
 
 
    ! allocate stuff
    ntype = this%ntype_of_particle
    safe_allocate(this%labels_particles(1:this%nparticle))
    safe_allocate(this%particletype(1:this%nparticle))
    safe_allocate(this%mass_particle(1:this%nparticle))
    safe_allocate(this%charge_particle(1:this%nparticle))
    safe_allocate(this%bosonfermion(1:this%nparticle))
    safe_allocate(this%nparticles_per_type(1:ntype))
    safe_allocate(this%particles_of_type(1:this%nparticle, 1:ntype))
 
    ! default all particles are electrons
    this%labels_particles = 'electron'
    this%particletype = 1
    this%mass_particle = m_one
    this%charge_particle = m_one
    this%bosonfermion = 1 ! set to fermion
 
    !%Variable DescribeParticlesModelmb
    !%Type block
    !%Section States::ModelMB
    !%Description
    !% Characterization of different modelmb particles in space%dim dimensional space.
    !%
    !% <tt>%DescribeParticlesModelmb
    !% <br>&nbsp;&nbsp; "proton"   | 1 | 1800. | 1. | fermion
    !% <br>&nbsp;&nbsp; "proton"   | 1 | 1800. | 1. | fermion
    !% <br>&nbsp;&nbsp; "electron" | 2 | 1.    | 1. | fermion
    !% <br>%</tt>
    !%
    !% would tell <tt>Octopus</tt> that there are presently 3 particles, called proton, proton,
    !% and electron, with types 1, 1, and 2, and corresponding masses and charges.
    !% All particles should be fermions, and this can be later enforced on the spatial
    !% part of the wavefunctions.
    !% The label and charge are presently only for informational purposes and
    !% are not checked or used in <tt>Octopus</tt>. The interaction has to take the
    !% actual charge into account.
    !%
    !%Option fermion 1
    !%  Particle is a fermion.
    !%Option boson 2
    !%  Particle is a boson.
    !%Option anyon 3
    !%  Particle is neither fermion nor boson.
    !%End
    if (parse_block(namespace, 'DescribeParticlesModelmb', blk) == 0) then
 
      call messages_experimental("Model many-body", namespace=namespace)
 
      ncols = parse_block_cols(blk, 0)
      if (ncols /= 5) then
        call messages_input_error(namespace, "DescribeParticlesModelmb")
      end if
      nline = parse_block_n(blk)
      if (nline /= this%nparticle) then
        call messages_input_error(namespace, "DescribeParticlesModelmb")
      end if
 
      do ipart = 1, this%nparticle
        call parse_block_string(blk, ipart - 1, 0, this%labels_particles(ipart))
        call parse_block_integer(blk, ipart - 1, 1, this%particletype(ipart))
        call parse_block_float(blk, ipart - 1, 2, this%mass_particle(ipart))
        call parse_block_float(blk, ipart - 1, 3, this%charge_particle(ipart))
        call parse_block_integer(blk, ipart - 1, 4, this%bosonfermion(ipart))
 
        write (message(1),'(a,a)') 'labels_particles = ', this%labels_particles(ipart)
        write (message(2),'(a,i6)') 'particletype = ', this%particletype(ipart)
        write (message(3),'(a,E20.10)') 'mass_particle = ', this%mass_particle(ipart)
        write (message(4),'(a,E20.10)') 'charge_particle = ', this%charge_particle(ipart)
        write (message(5),'(a,i6)') 'bosonfermion = ', this%bosonfermion(ipart)
        call messages_info(5, namespace=namespace)
      end do
      call parse_block_end(blk)
 
    end if
 
    this%nparticles_per_type = 0
    this%particles_of_type = 0
    do ipart = 1, this%nparticle
      this%nparticles_per_type(this%particletype(ipart)) = &
        this%nparticles_per_type(this%particletype(ipart)) + 1
      this%particles_of_type(this%nparticles_per_type(this%particletype(ipart)), &
        this%particletype(ipart)) = ipart
    end do
 
    this%max_particles_per_type = maxval(this%nparticles_per_type)
    itmp = this%max_particles_per_type
    jtmp = this%ntype_of_particle
    safe_allocate(this%exchange_symmetry(1:itmp, 1:itmp, 1:jtmp))
    this%exchange_symmetry = 0
 
    if (this%nparticle > 0) then
      ! FIXME: check why this is not initialized properly in the test, or why it is written out when not initialized
      safe_allocate(this%nspindown(1:this%ntype_of_particle, 1:nst))
      this%nspindown(:,:) = -1
      safe_allocate(this%iyoung(1:this%ntype_of_particle, 1:nst))
      this%iyoung(:,:) = -1
      safe_allocate(this%proj(1:nst))
      this%proj(:) = m_zero
    end if
 
    pop_sub(modelmb_particles_init)
 
  end subroutine modelmb_particles_init
 
 
  subroutine modelmb_particles_end (this)
    type(modelmb_particle_t),intent(inout) :: this
 
    push_sub(modelmb_particles_end)
 
    safe_deallocate_a(this%nspindown)
    safe_deallocate_a(this%iyoung)
    safe_deallocate_a(this%proj)
 
    safe_deallocate_a(this%labels_particles)
    safe_deallocate_a(this%particletype)
    safe_deallocate_a(this%mass_particle)
    safe_deallocate_a(this%charge_particle)
    safe_deallocate_a(this%nparticles_per_type)
    safe_deallocate_a(this%particles_of_type)
    safe_deallocate_a(this%exchange_symmetry)
    safe_deallocate_a(this%bosonfermion)
 
    safe_deallocate_a(this%labels_densities)
    safe_deallocate_a(this%particle_kept_densities)
 
    pop_sub(modelmb_particles_end)
  end subroutine modelmb_particles_end
 
  subroutine modelmb_particles_copy(modelmb_out, modelmb_in)
    type(modelmb_particle_t), intent(in)    :: modelmb_in
    type(modelmb_particle_t), intent(inout) :: modelmb_out
 
    push_sub(modelmb_particles_copy)
 
    call modelmb_particles_end(modelmb_out)
 
    modelmb_out%ndim = modelmb_in%ndim
    modelmb_out%ntype_of_particle = modelmb_in%ntype_of_particle
    modelmb_out%max_particles_per_type = modelmb_in%max_particles_per_type
    modelmb_out%nparticle = modelmb_in%nparticle
    modelmb_out%ndensities_to_calculate = modelmb_in%ndensities_to_calculate
 
    safe_allocate_source_a(modelmb_out%labels_particles,modelmb_in%labels_particles)
    safe_allocate_source_a(modelmb_out%particletype,modelmb_in%particletype)
    safe_allocate_source_a(modelmb_out%mass_particle,modelmb_in%mass_particle)
    safe_allocate_source_a(modelmb_out%charge_particle,modelmb_in%charge_particle)
    safe_allocate_source_a(modelmb_out%nparticles_per_type,modelmb_in%nparticles_per_type)
    safe_allocate_source_a(modelmb_out%particles_of_type,modelmb_in%particles_of_type)
    safe_allocate_source_a(modelmb_out%exchange_symmetry,modelmb_in%exchange_symmetry)
    safe_allocate_source_a(modelmb_out%bosonfermion,modelmb_in%bosonfermion)
 
    safe_allocate_source_a(modelmb_out%labels_densities,modelmb_in%labels_densities)
    safe_allocate_source_a(modelmb_out%particle_kept_densities,modelmb_in%particle_kept_densities)
 
    if (modelmb_in%nparticle > 0) then
      safe_allocate_source_a(modelmb_out%nspindown, modelmb_in%nspindown)
      safe_allocate_source_a(modelmb_out%iyoung, modelmb_in%iyoung)
      safe_allocate_source_a(modelmb_out%proj, modelmb_in%proj)
    end if
 
    pop_sub(modelmb_particles_copy)
 
  end subroutine modelmb_particles_copy
 
  subroutine modelmb_copy_masses(this, masses)
    class(modelmb_particle_t), intent(in)    :: this
    real(real64),             intent(inout) :: masses(:)
 
    real(real64), parameter :: tol_mass = 1.e-10_real64
    integer :: dimcounter,ipart
 
    push_sub(modelmb_copy_masses)
 
    ! copy masses to gr%der%masses
    dimcounter = 0
    do ipart = 1,this%nparticle
      if (abs(this%mass_particle(ipart)-1.0_real64) > tol_mass) then
        masses(dimcounter+1:dimcounter+this%ndim) = this%mass_particle(ipart)
      end if
      dimcounter = dimcounter+this%ndim
    end do
 
    pop_sub(modelmb_copy_masses)
  end subroutine modelmb_copy_masses
 
end module modelmb_particles_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: