LCOV - code coverage report
Current view: top level - src/78_effpot - m_lwf_primitive.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 75.2 % 234 176
Test Date: 2026-11-10 12:08:57 Functions: 77.8 % 9 7
Branches: 44.6 % 500 223

             Branch data     Line data    Source code
       1                 :             : !!****m*ABINIT/m_lwf_primitive_potential
       2                 :             : !! NAME
       3                 :             : !! m_lwf_primitive_potential
       4                 :             : !!
       5                 :             : !! FUNCTION
       6                 :             : !! This module contains the example lwf primitive potential. It is not the
       7                 :             : !! one really used. The purpose of this is to show how to extend multibinit
       8                 :             : !! to new degree of freedom
       9                 :             : !!
      10                 :             : !! Datatypes:
      11                 :             : !!  LWF_primitive_potential_t
      12                 :             : !!
      13                 :             : !! Subroutines:
      14                 :             : !!
      15                 :             : !! COPYRIGHT
      16                 :             : !! Copyright (C) 2001-2026 ABINIT group (hexu)
      17                 :             : !! This file is distributed under the terms of the
      18                 :             : !! GNU General Public License, see ~abinit/COPYING
      19                 :             : !! or http://www.gnu.org/copyleft/gpl.txt .
      20                 :             : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
      21                 :             : !!
      22                 :             : !! SOURCE
      23                 :             : 
      24                 :             : 
      25                 :             : #if defined HAVE_CONFIG_H
      26                 :             : #include "config.h"
      27                 :             : #endif
      28                 :             : #include "abi_common.h"
      29                 :             : 
      30                 :             : #define HAVE_NETCDF 1
      31                 :             : 
      32                 :             : module m_lwf_primitive_potential
      33                 :             :   use, intrinsic :: iso_c_binding
      34                 :             :   use defs_basis
      35                 :             :   use m_abicore
      36                 :             :   use m_errors
      37                 :             :   use m_nctk
      38                 :             : #if defined HAVE_NETCDF
      39                 :             :   use netcdf
      40                 :             : #endif
      41                 :             :   use m_xmpi
      42                 :             :   use m_mathfuncs, only: eigensh
      43                 :             :   use m_multibinit_dataset, only: multibinit_dtset_type
      44                 :             :   use m_multibinit_cell, only: mbcell_t, mbsupercell_t
      45                 :             :   use m_primitive_potential, only: primitive_potential_t
      46                 :             :   use m_abstract_potential, only: abstract_potential_t
      47                 :             :   use m_dynamic_array, only: int2d_array_type
      48                 :             :   use m_supercell_maker, only: supercell_maker_t
      49                 :             :   use m_spmat_ndcoo, only: ndcoo_mat_t
      50                 :             :   use m_lwf_potential, only: lwf_potential_t
      51                 :             :   implicit none
      52                 :             :   private
      53                 :             :   !!***
      54                 :             : 
      55                 :             :   !-------------------------------------------------------------------!
      56                 :             :   ! An harmonic potential which has only the IFC (with no dipole-dipole)
      57                 :             :   !
      58                 :             :   ! IFC is written in a coefficient matrix M(R, i, j)=val
      59                 :             :   !-------------------------------------------------------------------!
      60                 :             :   type, public, extends(primitive_potential_t):: lwf_primitive_potential_t
      61                 :             :      integer:: nlwf = 0, natom = 0, nR = 0  ! number of LWF
      62                 :             :      type(ndcoo_mat_t):: coeff  !  A N-dimensional COO matrix.
      63                 :             :      !  The indices are (ind_R, i, j). Note that xyz is included in i and j.
      64                 :             :      real(dp), allocatable:: lwf_masses(:)
      65                 :             : 
      66                 :             :      integer:: onebody_nterm = 0
      67                 :             : 
      68                 :             :      real(dp), allocatable:: lattice_coeffs(:,:, :) ! (natom*3, nlwf, nR)
      69                 :             : 
      70                 :             :      integer, allocatable:: onebody_i(:), onebody_order(:)
      71                 :             :      real(dp), allocatable:: onebody_val(:)
      72                 :             : 
      73                 :             :      type(ndcoo_mat_t):: coeff_twobody  ! Two body interaction parameters of higher order.
      74                 :             :      ! The indices are (indR, i, j, orderi, orderj), H(R, i, j)= val A_i^(orderi) A_j^(orderj)
      75                 :             :      integer, allocatable:: Rlist(:,:)  ! The list of R points (3, number of R-points)
      76                 :             : 
      77                 :             :      ! For testing the bounding.
      78                 :             :      !logical:: has_self_bound_term= .False.
      79                 :             :      !integer:: self_bound_order = 0
      80                 :             :      !real(dp):: self_bound_coeff = 0.0_dp
      81                 :             : 
      82                 :             :      !. Rlist(:, ind_R) is a R-vector.
      83                 :             :      real(dp):: ref_energy = 0.0                  ! reference energy
      84                 :             : 
      85                 :             :      logical:: as_lattice_anharmonic=.False.
      86                 :             : 
      87                 :             :    contains
      88                 :             :      procedure:: initialize
      89                 :             :      procedure:: finalize
      90                 :             :      procedure:: load_from_files   ! load potential from files listed in files file
      91                 :             :      procedure:: load_from_netcdf  ! load potential from a netcdf file
      92                 :             :      procedure:: fill_supercell     ! fill a supercell potential
      93                 :             :      procedure:: get_hamk          ! generate hamiltonian for one k point.
      94                 :             :      procedure:: get_eigen         ! eigen values and eigen vectors
      95                 :             :      !procedure:: add_self_bound_term  ! add a self bound term
      96                 :             :   end type lwf_primitive_potential_t
      97                 :             : 
      98                 :             : 
      99                 :             : contains
     100                 :             : 
     101                 :             :   !-------------------------------------------------------------------!
     102                 :             :   ! Initialize
     103                 :             :   !  Input:
     104                 :             :   !    nlwf: number of Wannier functions
     105                 :             :   !    primcell: the reference primitive cell.
     106                 :             :   !-------------------------------------------------------------------!
     107                 :           1 :   subroutine initialize(self, primcell)
     108                 :             :     class(lwf_primitive_potential_t), intent(inout):: self
     109                 :             :     type(mbcell_t), target, intent(inout):: primcell
     110                 :             :     !integer, intent(in):: nlwf
     111                 :           1 :     self%primcell=>primcell
     112                 :           1 :     self%label="lwf_primitive_potential"
     113                 :           1 :     self%has_spin=.False.
     114                 :           1 :     self%has_displacement=.False.
     115                 :           1 :     self%has_strain=.False.
     116                 :           1 :     self%has_lwf=.True.
     117                 :           1 :   end subroutine initialize
     118                 :             : 
     119                 :             : 
     120                 :             :   !-------------------------------------------------------------------!
     121                 :             :   ! Finalize
     122                 :             :   !-------------------------------------------------------------------!
     123                 :           1 :   subroutine finalize(self)
     124                 :             :     class(lwf_primitive_potential_t), intent(inout):: self
     125                 :           1 :     call self%coeff%finalize()
     126                 :           1 :     call self%coeff_twobody%finalize()
     127         [ +  - ]:           1 :     ABI_SFREE(self%Rlist)
     128                 :           1 :     nullify(self%primcell)
     129                 :           1 :     self%nlwf = 0
     130                 :           1 :     self%label="Destroyed lwf_primitive_potential"
     131                 :           1 :     call self%primitive_potential_t%finalize()
     132                 :             : 
     133         [ +  - ]:           1 :     ABI_SFREE(self%lattice_coeffs)
     134         [ +  - ]:           1 :     ABI_SFREE(self%lwf_masses)
     135                 :           1 :     self%natom = 0
     136                 :           1 :     self%nR = 0
     137                 :             : 
     138         [ +  - ]:           1 :     if (self%onebody_nterm /= 0) then
     139         [ +  - ]:           1 :        ABI_SFREE(self%onebody_i)
     140         [ +  - ]:           1 :        ABI_SFREE(self%onebody_val)
     141         [ +  - ]:           1 :        ABI_SFREE(self%onebody_order)
     142                 :             :     end if
     143                 :           1 :     self%onebody_nterm = 0
     144                 :           1 :   end subroutine finalize
     145                 :             : 
     146                 :             : 
     147                 :             :   !-------------------------------------------------------------------!
     148                 :             :   ! Load from files:
     149                 :             :   !> params: parameters
     150                 :             :   !> fnames: file names from files file
     151                 :             :   !-------------------------------------------------------------------!
     152                 :           1 :   subroutine load_from_files(self, params, fnames)
     153                 :             :     class(lwf_primitive_potential_t), intent(inout):: self
     154                 :             :     type(multibinit_dtset_type), intent(in):: params
     155                 :             :     character(len = fnlen), intent(in):: fnames(:)
     156                 :           1 :     call self%load_from_netcdf(fnames(1))
     157                 :             :     !call self%add_self_bound_term(params%lwf_self_bound_order, &
     158                 :             :     !     & params%lwf_self_bound_coeff)
     159                 :           1 :     self%as_lattice_anharmonic= (params%latt_lwf_anharmonic == 1)
     160                 :           1 :   end subroutine load_from_files
     161                 :             : 
     162                 :             :   !-------------------------------------------------------------------!
     163                 :             :   ! load potential from netcdf file
     164                 :             :   ! Note that the lattic part of the primitive cell is also loaded.
     165                 :             :   ! Input:
     166                 :             :   !  fname: filename
     167                 :             :   !-------------------------------------------------------------------!
     168                 :           1 :   subroutine load_from_netcdf(self, fname)
     169                 :             :     class(lwf_primitive_potential_t), intent(inout):: self
     170                 :             :     character(len = fnlen), intent(in):: fname
     171                 :             :     integer:: ncid, ierr
     172                 :             :     integer:: iR,  nR, nlwf, natom, twobody_nterm, onebody_nterm
     173                 :             :     real(dp):: cell(3, 3)
     174                 :           1 :     real(dp), allocatable:: ifc(:, :, :), xcart(:,:), masses(:)
     175                 :           1 :     real(dp), allocatable ::  twobody_val(:)
     176                 :           1 :     integer, allocatable ::   zion(:), twobody_iR(:), twobody_i(:), twobody_j(:), twobody_orderi(:), twobody_orderj(:)
     177                 :             :     integer:: varid, i, j
     178                 :             : #if defined HAVE_NETCDF
     179                 :           1 :     ierr = nf90_open(trim(fname), NF90_NOWRITE, ncid)
     180         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "Open netcdf file")
     181                 :             : 
     182                 :             :     ! read primcell info
     183                 :           1 :     ierr = nctk_get_dim(ncid, "wann_natom", natom)
     184         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "getting natom in lwf potential file")
     185                 :             : 
     186                 :           1 :     ierr = nctk_get_dim(ncid, "nR" , nR)
     187         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "getting nR in lwf potential file")
     188                 :           1 :     self%nR = nR
     189                 :             : 
     190                 :           1 :     ierr = nctk_get_dim(ncid, "wann_nwann", nlwf)
     191         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "getting wann_nwann in lwf potential file")
     192                 :             : 
     193                 :             :     ! TODO : add lattice
     194                 :             :     !call self%primcell%set_lattice()
     195                 :             : 
     196                 :           1 :     call self%primcell%set_lwf(nlwf)
     197                 :           1 :     self%nlwf = nlwf
     198         [ +  + ]:           4 :     call self%coeff%initialize([ nR, nlwf, nlwf])
     199         [ +  + ]:           6 :     call self%coeff_twobody%initialize([nR, nlwf, nlwf, -1, -1])
     200                 :             : 
     201   [ -  +  +  -  :           5 :     ABI_MALLOC(ifc, (nlwf, nlwf, nR))
          -  +  +  -  -  
          +  +  -  +  -  
          +  -  -  +  -  
                      + ]
     202   [ -  +  +  -  :           3 :     ABI_MALLOC(self%Rlist, (3, nR))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     203   [ -  +  +  -  :           3 :     ABI_MALLOC(xcart, (3, natom))
          +  -  +  -  -  
                +  -  + ]
     204   [ -  +  +  -  :           3 :     ABI_MALLOC(zion, (natom))
          +  -  +  -  -  
                +  -  + ]
     205   [ -  +  +  -  :           3 :     ABI_MALLOC(masses, (natom))
          +  -  +  -  -  
                +  -  + ]
     206                 :             : 
     207   [ -  +  +  -  :           3 :     ABI_MALLOC(self%lwf_masses, (nlwf))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     208                 :             : 
     209                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_lwf_masses", varid)
     210         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "lwf_masses")
     211                 :           1 :     ierr = nf90_get_var(ncid, varid, self%lwf_masses)
     212         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "lwf_masses")
     213         [ +  + ]:           3 :     self%lwf_masses(:) = self%lwf_masses(:) * amu_emass
     214   [ +  -  -  +  :           4 :     self%primcell%lwf%lwf_masses = self%lwf_masses
          -  -  -  -  -  
                -  +  + ]
     215                 :             : 
     216                 :           1 :     ierr = nf90_inq_dimid(ncid, "wann_onebody_nterm", onebody_nterm)
     217         [ -  + ]:           1 :     if (ierr /= nf90_noerr) then
     218                 :           0 :        onebody_nterm = 0
     219                 :             :     else
     220                 :           1 :        ierr = nctk_get_dim(ncid, "wann_onebody_nterm", onebody_nterm)
     221         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "getting wann_onebody_nterm in lwf potential file")
     222                 :             :     end if
     223                 :           1 :     self%onebody_nterm = onebody_nterm
     224                 :             : 
     225                 :             : 
     226                 :           1 :     ierr = nf90_inq_dimid(ncid, "wann_twobody_nterm", twobody_nterm)
     227         [ +  - ]:           1 :     if (ierr /= nf90_noerr) then
     228                 :           1 :        twobody_nterm = 0
     229                 :             :     else
     230                 :           0 :        ierr = nctk_get_dim(ncid, "wann_twobody_nterm", twobody_nterm)
     231         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "getting wann_twobody_nterm in lwf potential file")
     232                 :             :     end if
     233                 :             : 
     234                 :           1 :     self%natom = natom
     235                 :             : 
     236                 :             : 
     237                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_cell", varid)
     238         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_cell")
     239                 :           1 :     ierr = nf90_get_var(ncid, varid, cell)
     240         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_cell")
     241   [ +  +  +  + ]:          13 :     cell(:,:)=cell(:,:)/ Bohr_Ang
     242                 :             : 
     243                 :             : 
     244                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_atomic_xcart", varid)
     245         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_atomic_xcart")
     246                 :           1 :     ierr = nf90_get_var(ncid, varid, xcart)
     247         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_atomic_xcart")
     248   [ +  +  +  + ]:          25 :     xcart(:,:)=xcart(:,:)/ Bohr_Ang
     249                 :             : 
     250                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_atomic_numbers", varid)
     251         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_atomic_numbers")
     252                 :           1 :     ierr = nf90_get_var(ncid, varid, zion)
     253         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_atomic_numbers")
     254                 :             : 
     255                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_atomic_masses", varid)
     256         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_atomic_masses")
     257                 :           1 :     ierr = nf90_get_var(ncid, varid, masses)
     258         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_atomic_masses")
     259                 :             : 
     260                 :           1 :     call self%primcell%set_lattice(natom, cell, xcart, masses, zion)
     261                 :             : 
     262                 :             : 
     263   [ -  +  +  -  :           5 :     ABI_MALLOC(self%lattice_coeffs, (nlwf, natom*3, nR))
          +  -  -  +  +  
          -  -  +  +  -  
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     264                 :             : 
     265                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_wannier_function_real", varid)
     266         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_wannier_function_real")
     267                 :           1 :     ierr = nf90_get_var(ncid, varid, self%lattice_coeffs)
     268         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_wannier_function_real")
     269                 :             : 
     270                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_Rlist", varid)
     271         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_Rlist")
     272                 :           1 :     ierr = nf90_get_var(ncid, varid, self%Rlist)
     273         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_Rlist")
     274                 :             : 
     275                 :           1 :     ierr = nf90_inq_varid(ncid, "wann_HamR_real", varid)
     276         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_HamR_real")
     277                 :           1 :     ierr = nf90_get_var(ncid, varid, ifc)
     278         [ -  + ]:           1 :     NCF_CHECK_MSG(ierr, "wann_HamR_real")
     279                 :             : 
     280                 :             : 
     281   [ +  +  +  +  :         316 :     ifc(:,:,:) = ifc(:,:,:) * eV_Ha * (Bohr_Ang*Bohr_Ang)
                   +  + ]
     282                 :             : 
     283         [ +  + ]:          46 :     do iR = 1, nR
     284         [ +  + ]:         136 :        do i = 1, nlwf
     285         [ +  + ]:         315 :           do j = 1, nlwf
     286         [ +  + ]:         270 :              if (abs(ifc(j, i, iR))>1e-4) then
     287                 :             :                 ! NOTE: in fortran the order of index in reversed when reading netcdf array.
     288         [ +  + ]:         528 :                 call self%coeff%add_entry([iR, i, j], ifc(j, i, iR))
     289                 :             :              end if
     290                 :             :           end do
     291                 :             :        end do
     292                 :             :     end do
     293                 :             : 
     294                 :             : 
     295         [ +  - ]:           1 :     if(self%onebody_nterm /= 0) then
     296   [ -  +  +  -  :           3 :        ABI_MALLOC(self%onebody_i, (self%onebody_nterm))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     297   [ -  +  +  -  :           2 :        ABI_MALLOC(self%onebody_order, (self%onebody_nterm))
          -  +  -  +  -  
                      + ]
     298   [ -  +  +  -  :           3 :        ABI_MALLOC(self%onebody_val, (self%onebody_nterm))
          +  -  +  -  -  
             +  -  +  -  
                      + ]
     299                 :             : 
     300                 :           1 :        ierr = nf90_inq_varid(ncid, "wann_onebody_i", varid)
     301         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "wann_onebody_i")
     302                 :           1 :        ierr = nf90_get_var(ncid, varid, self%onebody_i)
     303         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "wann_onebody_i")
     304                 :             : 
     305         [ +  + ]:           7 :        self%onebody_i = self%onebody_i+1
     306                 :             : 
     307                 :           1 :        ierr = nf90_inq_varid(ncid, "wann_onebody_val", varid)
     308         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "wann_onebody_val")
     309                 :           1 :        ierr = nf90_get_var(ncid, varid, self%onebody_val)
     310         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "wann_onebody_val")
     311                 :             : 
     312                 :           1 :        ierr = nf90_inq_varid(ncid, "wann_onebody_order", varid)
     313         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "wann_onebody_order")
     314                 :           1 :        ierr = nf90_get_var(ncid, varid, self%onebody_order)
     315         [ -  + ]:           1 :        NCF_CHECK_MSG(ierr, "wann_onebody_order")
     316                 :             : 
     317         [ +  + ]:           7 :        do i = 1, self%onebody_nterm
     318                 :           7 :           self%onebody_val(i)=self%onebody_val(i) * eV_Ha * (Bohr_Ang**self%onebody_order(i))
     319                 :             :        end do
     320                 :             :     end if
     321                 :             : 
     322         [ -  + ]:           1 :     if(twobody_nterm /= 0) then
     323   [ #  #  #  #  :           0 :        ABI_MALLOC(twobody_iR, (twobody_nterm))
          #  #  #  #  #  
                #  #  # ]
     324   [ #  #  #  #  :           0 :        ABI_MALLOC(twobody_i, (twobody_nterm))
             #  #  #  # ]
     325   [ #  #  #  #  :           0 :        ABI_MALLOC(twobody_j, (twobody_nterm))
             #  #  #  # ]
     326   [ #  #  #  #  :           0 :        ABI_MALLOC(twobody_orderi, (twobody_nterm))
             #  #  #  # ]
     327   [ #  #  #  #  :           0 :        ABI_MALLOC(twobody_orderj, (twobody_nterm))
             #  #  #  # ]
     328   [ #  #  #  #  :           0 :        ABI_MALLOC(twobody_val, (twobody_nterm))
          #  #  #  #  #  
                #  #  # ]
     329                 :             : 
     330                 :           0 :        ierr = nf90_inq_varid(ncid, "wann_twobody_iR", varid)
     331         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_iR")
     332                 :           0 :        ierr = nf90_get_var(ncid, varid, twobody_iR)
     333         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_iR")
     334                 :             : 
     335                 :           0 :        ierr = nf90_inq_varid(ncid, "wann_twobody_i", varid)
     336         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_i")
     337                 :           0 :        ierr = nf90_get_var(ncid, varid, twobody_i)
     338         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_i")
     339                 :             : 
     340                 :           0 :        ierr = nf90_inq_varid(ncid, "wann_twobody_j", varid)
     341         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_j")
     342                 :           0 :        ierr = nf90_get_var(ncid, varid, twobody_j)
     343         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_j")
     344                 :             : 
     345                 :           0 :        ierr = nf90_inq_varid(ncid, "wann_twobody_orderi", varid)
     346         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_orderi")
     347                 :           0 :        ierr = nf90_get_var(ncid, varid, twobody_orderi)
     348         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_i")
     349                 :             : 
     350                 :           0 :        ierr = nf90_inq_varid(ncid, "wann_twobody_orderj", varid)
     351         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_orderj")
     352                 :           0 :        ierr = nf90_get_var(ncid, varid, twobody_orderj)
     353         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_orderj")
     354                 :             : 
     355                 :           0 :        ierr = nf90_inq_varid(ncid, "wann_twobody_val", varid)
     356         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_val")
     357                 :           0 :        ierr = nf90_get_var(ncid, varid, twobody_val)
     358         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "wann_twobody_val")
     359                 :             : 
     360                 :           0 :        ierr = nf90_close(ncid)
     361         [ #  # ]:           0 :        NCF_CHECK_MSG(ierr, "Close netcdf file")
     362                 :             : 
     363         [ #  # ]:           0 :        do i = 1, twobody_nterm
     364                 :             :         twobody_val(i)=twobody_val(i) * eV_Ha *  &
     365                 :           0 :             & (Bohr_Ang**(twobody_orderi(i)+ twobody_orderj(i)))
     366                 :             :           call self%coeff_twobody%add_entry([twobody_iR(i), twobody_i(i), &
     367         [ #  # ]:           0 :               & twobody_j(i), twobody_orderi(i), twobody_orderj(i)], twobody_val(i))
     368                 :             :        end do
     369                 :             :     end if
     370                 :             : 
     371         [ +  - ]:           1 :     ABI_SFREE(xcart)
     372         [ +  - ]:           1 :     ABI_SFREE(zion)
     373         [ +  - ]:           1 :     ABI_SFREE(masses)
     374         [ +  - ]:           1 :     ABI_SFREE(ifc)
     375                 :             : 
     376         [ -  + ]:           1 :     if(twobody_nterm /= 0) then
     377         [ #  # ]:           0 :        ABI_SFREE(twobody_i)
     378         [ #  # ]:           0 :        ABI_SFREE(twobody_j)
     379         [ #  # ]:           0 :        ABI_SFREE(twobody_orderi)
     380         [ #  # ]:           0 :        ABI_SFREE(twobody_orderj)
     381         [ #  # ]:           0 :        ABI_SFREE(twobody_val)
     382                 :             :     endif
     383                 :             : #else
     384                 :             :     NETCDF_NOTENABLED_ERROR()
     385                 :             : #endif
     386                 :             : 
     387   [ -  +  -  +  :           1 :   end subroutine load_from_netcdf
          -  +  -  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
                      + ]
     388                 :             : 
     389                 :             :   !-------------------------------------------------------------------!
     390                 :             :   !Fill supercell
     391                 :             :   ! Inputs:
     392                 :             :   !  scmaker : supercell_maker_t
     393                 :             :   ! Output:
     394                 :             :   !  scpot: a class pointer to an ABSTRACT potential.
     395                 :             :   !         Not a  potential because this is inherited from
     396                 :             :   !         an abstract_primitive_potential_t, which doesn't know
     397                 :             :   !         the type of the supercell potential.
     398                 :             :   !
     399                 :             :   !-------------------------------------------------------------------!
     400                 :           1 :   subroutine fill_supercell(self, scmaker, params, scpot, supercell)
     401                 :             :     use m_spmat_convert, only: COO_to_dense
     402                 :             : 
     403                 :             :     class(lwf_primitive_potential_t), intent(inout):: self
     404                 :             :     type(supercell_maker_t),                        intent(inout):: scmaker
     405                 :             :     type(multibinit_dtset_type),                    intent(inout):: params
     406                 :             :     class(abstract_potential_t), pointer,           intent(inout):: scpot
     407                 :             :     type(mbsupercell_t), target:: supercell
     408                 :             : 
     409                 :             :     integer:: nlwf, sc_nlwf
     410                 :             :     integer:: inz, iR, R(3), i, j, icell
     411                 :           1 :     integer, allocatable:: ilist_sc(:), jlist_sc(:), Rlist_sc(:,:)
     412                 :             :     real(dp):: val
     413                 :             : 
     414                 :           1 :     ABI_UNUSED_A(params)
     415                 :             : 
     416                 :             : 
     417                 :           1 :     nlwf = self%nlwf
     418                 :           1 :     sc_nlwf = nlwf*scmaker%ncells
     419                 :             : 
     420                 :             :     ! Step 1: allocate the scpot as a corresponding supercell potential
     421         [ -  + ]:           1 :     ABI_MALLOC_TYPE_SCALAR(lwf_potential_t, scpot)
     422                 :             :     ! Fortran does not know the functions specific to the derived class pointer.
     423                 :             :     ! Only the ones inheritated from abstract class,
     424                 :             :     ! unless select type is used:
     425                 :             :     select type(scpot)
     426                 :             :     type is (lwf_potential_t)
     427                 :           1 :        call scpot%initialize(sc_nlwf)
     428                 :           1 :        call scpot%set_supercell(supercell)
     429                 :             :        ! IFC is an COO_mat_t, which has the index of R1, R2, R3, i, j and the value of val
     430                 :             :        ! list of index R: coeff%ind%data(1, 1:coeff%nnz)
     431                 :             :        ! list of i: coeff%ind%data(2, 1:coeff%nnz)
     432                 :             :        ! list of j: coeff%ind%data(3, 1:coeff%nnz)
     433                 :             :        ! IFC: (ind) = val
     434         [ -  + ]:           1 :        if(self%as_lattice_anharmonic) then
     435                 :           0 :           call scpot%use_as_lattice_anharmonic()
     436                 :             :        else
     437                 :             :           ! harmonic terms
     438         [ +  + ]:         133 :           do inz = 1, self%coeff%nnz
     439                 :             :              ! For each non-zero entry in the coeff matrix
     440                 :             :              ! get the R, i, and j, val
     441                 :         132 :              iR = self%coeff%ind%data(1, inz)
     442         [ +  + ]:         528 :              R(:) = self%Rlist(:, iR)
     443                 :         132 :              i = self%coeff%ind%data(2, inz)
     444                 :         132 :              j = self%coeff%ind%data(3, inz)
     445                 :         132 :              val = self%coeff%val%data(inz)
     446                 :             :              ! translate i to i in supercell.
     447                 :             :              ! No need to allocate, it is done by trans_i . but remember to deallocate!
     448                 :             :              ! nbasis is the number in one primitive cell.
     449                 :             :              ! e.g. there are 3*natom possible i (3: x, y, z) in each primitive cell.
     450                 :         132 :              call scmaker%trans_i(nbasis = self%nlwf, i = i, i_sc = ilist_sc )
     451                 :             :              ! translate j, Rj to supercell.
     452                 :         132 :              call scmaker%trans_j_and_Rj(nbasis = self%nlwf, j = j, Rj = R, j_sc = jlist_sc, Rj_sc = Rlist_sc)
     453                 :             :              ! values are repeated in cells
     454         [ +  + ]:       67716 :              do icell = 1, scmaker%ncells
     455                 :       67716 :                 call scpot%add_term(ilist_sc(icell), jlist_sc(icell), val )
     456                 :             :              end do
     457         [ +  - ]:         132 :              ABI_SFREE(ilist_sc)
     458         [ +  - ]:         132 :              ABI_SFREE(jlist_sc)
     459         [ +  - ]:         133 :              ABI_SFREE(Rlist_sc)
     460                 :             :           end do
     461                 :             :           end if
     462                 :             : 
     463                 :             :        ! coefficients of atomic displacements in supercell
     464                 :           1 :        call scpot%supercell%lwf%lwf_latt_coeffs%initialize(self%nlwf*scmaker%ncells, self%natom*3*scmaker%ncells )
     465         [ +  + ]:           3 :        do i = 1, self%nlwf
     466                 :           2 :           call scmaker%trans_i(nbasis = nlwf, i = i, i_sc = ilist_sc)
     467                 :             : 
     468         [ +  + ]:          92 :           do iR = 1, self%nR
     469         [ +  + ]:         360 :              R(:) = self%Rlist(:, iR)
     470         [ +  + ]:        1712 :              do j = 1, self%natom*3
     471                 :        1620 :                 val = self%lattice_coeffs(i, j, iR)  ! ilwf, iatom3, iR
     472         [ +  + ]:        1710 :                 if (abs(val) > 1e-4) then
     473                 :        1348 :                    call scmaker%trans_j_and_Rj(nbasis = self%natom*3, j = j, Rj = R, j_sc = jlist_sc, Rj_sc = Rlist_sc)
     474         [ +  + ]:      691524 :                    do icell = 1, scmaker%ncells
     475                 :             :                       !call scpot%supercell%lwf%lwf_latt_coeffs(ilist_sc(icell))%push(jlist_sc(icell), val)
     476         [ +  + ]:     2071876 :                       call scpot%supercell%lwf%lwf_latt_coeffs%coeffs%add_entry([jlist_sc(icell), ilist_sc(icell)], val)
     477                 :             :                    end do
     478         [ +  - ]:        1348 :                    ABI_SFREE(jlist_sc)
     479         [ +  - ]:        1348 :                    ABI_SFREE(Rlist_sc)
     480                 :             :                 end if
     481                 :             :              end do
     482                 :             :           end do
     483         [ +  - ]:           3 :           ABI_SFREE(ilist_sc)
     484                 :             :        end do
     485                 :             : 
     486                 :             :        ! anharmonic terms
     487                 :             :        !call scpot%set_ref_energy(self%ref_energy*scmaker%ncells)
     488                 :             :        !if (self%has_self_bound_term) then
     489                 :             :        !   call scpot%add_self_bound_term(self%self_bound_order, self%self_bound_coeff)
     490                 :             :        !end if
     491                 :             : 
     492         [ +  - ]:           2 :        if (self%onebody_nterm /= 0) then
     493         [ +  + ]:           7 :           do i = 1, self%onebody_nterm
     494                 :           6 :              call scmaker%trans_i(nbasis = self%nlwf, i = self%onebody_i(i), i_sc = ilist_sc)
     495         [ +  + ]:        3078 :              do icell = 1, scmaker%ncells
     496                 :        3078 :                 call scpot%add_onebody_term(i = ilist_sc(icell), order = self%onebody_order(i), val = self%onebody_val(i))
     497                 :             :              end do
     498         [ +  - ]:           7 :              ABI_SFREE(ilist_sc)
     499                 :             :           end do
     500                 :             :        endif
     501                 :             :        ! Test the phonon energy
     502                 :             :        !call COO_to_dense(scpot%coeff, real_sc_evecs)
     503                 :             :        !sc_evecs(:,:) = real_sc_evecs(:,:)
     504                 :             :        !call eigensh(sc_evals, sc_evecs)
     505                 :             : 
     506                 :             :     end select
     507                 :             : 
     508                 :             : 
     509   [ -  +  -  +  :           1 :   end subroutine fill_supercell
                   -  + ]
     510                 :             : 
     511                 :             :   !-------------------------------------------------------------------!
     512                 :             :   ! calculate hamiltonian at k point (or more precisely, q point)
     513                 :             :   ! H_ij(k) = sum_R H_ij(R) exp( i2pi k.dot.R)
     514                 :             :   !-------------------------------------------------------------------!
     515         [ #  # ]:           0 :   subroutine get_hamk(self, kpoint, hamk)
     516                 :             :     class(lwf_primitive_potential_t), intent(in):: self
     517                 :             :     real(dp), intent(in):: kpoint(3)
     518                 :             :     complex(dp), intent(inout):: hamk(:,:)
     519                 :             :     integer:: inz, iR, R(3), i, j
     520                 :             :     real(dp):: val
     521   [ #  #  #  # ]:           0 :     hamk(:,:) = cmplx(0.0, 0.0)
     522         [ #  # ]:           0 :     do inz = 1, self%coeff%nnz
     523                 :           0 :        iR = self%coeff%ind%data(1, inz)
     524         [ #  # ]:           0 :        R(:) = self%Rlist(:, iR)
     525                 :           0 :        i = self%coeff%ind%data(2, inz)
     526                 :           0 :        j = self%coeff%ind%data(3, inz)
     527                 :           0 :        val = self%coeff%val%data(inz)
     528         [ #  # ]:           0 :        hamk(i, j) =hamk(i, j) + val*exp(cmplx(0.0, 2.0) *pi*dot_product(kpoint, R))
     529                 :             :     end do
     530                 :           0 :   end subroutine get_hamk
     531                 :             : 
     532                 :             :   !-------------------------------------------------------------------!
     533                 :             :   !calculate eigenvalue and eigen vector
     534                 :             :   !-------------------------------------------------------------------!
     535   [ #  #  #  # ]:           0 :   subroutine get_eigen(self, kpoint, evals, evecs)
     536                 :             :     class(lwf_primitive_potential_t), intent(in):: self
     537                 :             :     real(dp), intent(in):: kpoint(3)
     538                 :             :     real(dp), intent(inout):: evals(:)
     539                 :             :     complex(dp), intent(inout):: evecs(:,:)
     540                 :           0 :     call self%get_hamk(kpoint, evecs)
     541                 :             :     ! The evecs array is reused both as the matrix and eigenvectors.
     542                 :           0 :     call eigensh(evals, evecs)
     543                 :           0 :   end subroutine get_eigen
     544                 :             : 
     545                 :             :   !subroutine add_self_bound_term(self, order, coeff)
     546                 :             :   !  class(lwf_primitive_potential_t), intent(inout):: self
     547                 :             :   !  integer, intent(in):: order
     548                 :             :   !  real(dp), intent(in):: coeff
     549                 :             :   !  if (order /= 0) then
     550                 :             :   !     self%has_self_bound_term=.True.
     551                 :             :   !     self%self_bound_order = order
     552                 :             :   !     self%self_bound_coeff = coeff
     553                 :             :   !  end if
     554                 :             : 
     555                 :             :   !end subroutine add_self_bound_term
     556                 :             : 
     557                 :             : 
     558   [ +  +  -  +  :           4 : end module m_lwf_primitive_potential
          +  +  -  +  -  
          +  -  +  -  +  
          -  +  +  +  -  
          +  -  +  -  +  
          -  +  -  +  -  
          +  -  +  -  +  
             -  +  -  + ]
        

Generated by: LCOV version 2.3-1