TABLE OF CONTENTS


ABINIT/m_dens [ Modules ]

[ Top ] [ Modules ]

NAME

  m_dens

FUNCTION

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (MT,ILuk,MVer,EB,SPr)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_dens
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_xmpi
32  use m_splines
33 
34  use defs_abitypes, only : MPI_type
35  use m_time,        only : timab
36  use m_numeric_tools, only : wrap2_zero_one
37  use m_io_tools,    only : open_file
38  use m_geometry,    only : xcart2xred, metric
39  use m_mpinfo,      only : ptabs_fourdp
40 
41  implicit none
42 
43  private
44 
45  public :: dens_hirsh      ! Compute the Hirshfeld charges
46  public :: mag_constr      ! Compute the potential corresponding to constrained magnetic moments.
47  public :: mag_constr_e    ! Compute the energy corresponding to constrained magnetic moments.
48  public :: calcdensph      ! Compute and print integral of total density inside spheres around atoms.

ABINIT/printmagvtk [ Functions ]

[ Top ] [ Functions ]

NAME

  printmagvtk

FUNCTION

  Auxiliary routine for printing out magnetization density in VTK format.
  Output file name is DEN.vtk

INPUTS

  mpi_enreg = information about adopted parallelization strategy
  nspden    = number of components of density matrix (possible values re 1,2, or 4)
              nspden:   1 -> rho
              nspden:   2 -> rho_up,rho_dwn
              nspden:   4 -> rho,mx,my,mz
  nfft      = number of fft points per FFT processor
  ngfft     = full information about FFT mesh
  rhor      = density array stored in the memory of current FFT processor
  rprimd    = array of lattice vectors

OUTPUT

SIDE EFFECTS

NOTES

  At the moment this routine is mainly used for development and debugging
  of gs and dfpt calculations with non-collinear spins. If needed, can be used
  to print final density in vtk format.
  IMPORTANT: implementation is thoroughly checked only for npspinor = 1,
             for other case might need to change the part gathering
             the FFT mesh info

PARENTS

CHILDREN

      xmpi_gather

SOURCE

1436 subroutine printmagvtk(mpi_enreg,cplex,nspden,nfft,ngfft,rhor,rprimd,fname)
1437 
1438 
1439 !This section has been created automatically by the script Abilint (TD).
1440 !Do not modify the following lines by hand.
1441 #undef ABI_FUNC
1442 #define ABI_FUNC 'printmagvtk'
1443 !End of the abilint section
1444 
1445  implicit none
1446 
1447 !Arguments ------------------------------------
1448 !scalars
1449  type(MPI_type),intent(in)   :: mpi_enreg
1450  integer,intent(in)          :: nfft,nspden,cplex
1451 !arrays
1452  integer,intent(in)          :: ngfft(18)
1453  real(dp),intent(in)         :: rhor(cplex*nfft,nspden),rprimd(3,3)
1454  character(len=*),intent(in) :: fname
1455 
1456 !Local variables-------------------------------
1457 !scalars
1458  integer :: denvtk,denxyz,denxyz_im,nfields
1459  integer :: nx,ny,nz,nfft_tot
1460  integer :: ii,jj,kk,ind,ispden
1461  integer :: mpi_comm,mpi_head,mpi_rank,ierr
1462  real    :: rx,ry,rz
1463  integer :: nproc_fft,ir
1464  character(len=500) :: msg
1465  character(len=10)  :: outformat
1466  character(len=50)   :: fname_vtk
1467  character(len=50)   :: fname_xyz
1468  character(len=50)   :: fname_xyz_re
1469  character(len=50)   :: fname_xyz_im
1470 !arrays
1471  real(dp),allocatable :: rhorfull(:,:)
1472 
1473 
1474 ! *************************************************************************
1475 
1476  DBG_ENTER("COLL")
1477 
1478 ! if (option/=1 .and. option/=2 ) then
1479 !   write(msg,'(3a,i0)')&
1480 !&   'The argument option should be 1 or 2,',ch10,&
1481 !&   'however, option=',option
1482 !   MSG_BUG(msg)
1483 ! end if
1484 !
1485 ! if (sizein<1) then
1486 !   write(msg,'(3a,i0)')&
1487 !&   'The argument sizein should be a positive number,',ch10,&
1488 !&   'however, sizein=',sizein
1489 !   MSG_ERROR(msg)
1490 ! end if
1491 
1492  DBG_EXIT("COLL")
1493 
1494  fname_vtk=adjustl(adjustr(fname)//".vtk")
1495  fname_xyz=adjustl(adjustr(fname)//".xyz")
1496  fname_xyz_re=adjustl(adjustr(fname)//"_re.xyz")
1497  fname_xyz_im=adjustl(adjustr(fname)//"_im.xyz")
1498  !write(std_out,*) ' Writing out .vtk file: ',fname_vtk
1499  !write(std_out,*) ' Writing out .xyz file: ',fname_xyz
1500 
1501   !if 1 or two component density then write out either 1 or 2 scalar density fields
1502   !if 4, then write one scalar field (density) and one vector field (magnetization density)
1503  if(nspden/=4)then
1504    nfields=nspden
1505  else
1506    nfields=2
1507  end if
1508 
1509  nfields=nfields*cplex
1510 
1511   ! FFT mesh specifications: full grid
1512  nx=ngfft(1)        ! number of points along 1st lattice vector
1513  ny=ngfft(2)        ! number of points along 2nd lattice vector
1514  nz=ngfft(3)        ! number of points along 3rd lattice vector
1515  nfft_tot=nx*ny*nz  ! total number of fft mesh points (can be different from nfft in case of distributed memory of nproc_fft processors)
1516 
1517 
1518   ! Gather information about memory distribution
1519  mpi_head=0
1520  mpi_comm = mpi_enreg%comm_fft
1521  mpi_rank = xmpi_comm_rank(mpi_comm)
1522  nproc_fft=ngfft(10)
1523 
1524   ! Create array to host full FFT mesh
1525  if(mpi_rank==mpi_head)then
1526    ABI_ALLOCATE(rhorfull,(cplex*nfft_tot,nspden))
1527  end if
1528 
1529   ! Fill in the full mesh
1530  if(nproc_fft==1)then
1531    rhorfull=rhor
1532  else
1533    do ir=1,nspden
1534      call xmpi_gather(rhor(:,ir),cplex*nfft,rhorfull(:,ir),cplex*nfft,mpi_head,mpi_comm,ierr)
1535    end do
1536  end if
1537 
1538  if(mpi_rank==mpi_head)then
1539 
1540     ! Open the output vtk file
1541    if (open_file(fname_vtk,msg,newunit=denvtk,status='replace',form='formatted') /=0) then
1542      MSG_WARNING(msg)
1543      RETURN
1544    end if
1545 
1546    if(cplex==1) then
1547      if (open_file(fname_xyz,msg,newunit=denxyz,status='replace',form='formatted') /=0) then
1548        MSG_WARNING(msg)
1549        RETURN
1550      end if
1551    else if (cplex==2) then
1552      if (open_file(fname_xyz_re,msg,newunit=denxyz,status='replace',form='formatted') /=0) then
1553        MSG_WARNING(msg)
1554        RETURN
1555      end if
1556      if (open_file(fname_xyz_im,msg,newunit=denxyz_im,status='replace',form='formatted') /=0) then
1557        MSG_WARNING(msg)
1558        RETURN
1559      end if
1560    end if
1561 
1562     ! Write the header of the output vtk file
1563    write(denvtk,"(a)") '# vtk DataFile Version 2.0'
1564    write(denvtk,"(a)") 'Electron density components'
1565    write(denvtk,"(a)") 'ASCII'
1566    write(denvtk,"(a)") 'DATASET STRUCTURED_GRID'
1567    write(denvtk,"(a,3i6)") 'DIMENSIONS ', nx,ny,nz
1568    write(denvtk,"(a,i18,a)") 'POINTS ',nfft_tot,' double'
1569 
1570    if (nspden==1) then
1571      outformat="(4e16.8)"
1572    else if (nspden==2) then
1573      outformat="(5e16.8)"
1574    else
1575      outformat="(7e16.8)"
1576    end if
1577 
1578     ! Write out information about grid points
1579    do kk=0,nz-1
1580      do jj=0,ny-1
1581        do ii=0,nx-1
1582 
1583          rx=(dble(ii)/nx)*rprimd(1,1)+(dble(jj)/ny)*rprimd(1,2)+(dble(kk)/nz)*rprimd(1,3)
1584          ry=(dble(ii)/nx)*rprimd(2,1)+(dble(jj)/ny)*rprimd(2,2)+(dble(kk)/nz)*rprimd(2,3)
1585          rz=(dble(ii)/nx)*rprimd(3,1)+(dble(jj)/ny)*rprimd(3,2)+(dble(kk)/nz)*rprimd(3,3)
1586          write(denvtk,'(3f16.8)') rx,ry,rz  !coordinates of the grid point
1587          ind=1+ii+nx*(jj+ny*kk)
1588          if (cplex==1) then
1589            write(denxyz,outformat) rx,ry,rz,(rhorfull(ind,ispden),ispden=1,nspden)
1590          else
1591            write(denxyz,outformat)    rx,ry,rz,(rhorfull(2*ind-1,ispden),ispden=1,nspden)
1592            write(denxyz_im,outformat) rx,ry,rz,(rhorfull(2*ind  ,ispden),ispden=1,nspden)
1593          end if
1594        end do
1595      end do
1596    end do
1597 
1598    if(cplex==1) then
1599      close(denxyz)
1600    else
1601      close(denxyz)
1602      close(denxyz_im)
1603    end if
1604 
1605     ! Write out information about field defined on the FFT mesh
1606    write(denvtk,"(a,i18)") 'POINT_DATA ',nfft_tot
1607    write(denvtk,"(a,i6)")  'FIELD Densities ',nfields
1608 
1609 
1610     ! Write out different fields depending on the number of density matrix components
1611    if(nspden==1)then
1612 
1613       !single component, so just write out the density
1614      if(cplex==1) then
1615        write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
1616        do kk=0,nz-1
1617          do jj=0,ny-1
1618            do ii=0,nx-1
1619              ind=1+ii+nx*(jj+ny*kk)
1620              write(denvtk,'(f16.8)') rhorfull(ind,1)
1621            end do
1622          end do
1623        end do
1624      else
1625        write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
1626        do kk=0,nz-1
1627          do jj=0,ny-1
1628            do ii=0,nx-1
1629              ind=2*(1+ii+nx*(jj+ny*kk))-1
1630              write(denvtk,'(f16.8)') rhorfull(ind,1)
1631            end do
1632          end do
1633        end do
1634        write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double'
1635        do kk=0,nz-1
1636          do jj=0,ny-1
1637            do ii=0,nx-1
1638              ind=2*(1+ii+nx*(jj+ny*kk))
1639              write(denvtk,'(f16.8)') rhorfull(ind,1)
1640            end do
1641          end do
1642        end do
1643      end if
1644 
1645    else if(nspden==2)then
1646 
1647       !two component, write the density for spin_up and spin_down channels
1648      if(cplex==1) then
1649 
1650        write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
1651        do kk=0,nz-1
1652          do jj=0,ny-1
1653            do ii=0,nx-1
1654              ind=1+ii+nx*(jj+ny*kk)
1655              write(denvtk,'(f16.8)') rhorfull(ind,1)
1656            end do
1657          end do
1658        end do
1659        write(denvtk,"(a,i18,a)") 'mag 1 ',nfft_tot,' double'
1660        do kk=0,nz-1
1661          do jj=0,ny-1
1662            do ii=0,nx-1
1663              ind=1+ii+nx*(jj+ny*kk)
1664              write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
1665            end do
1666          end do
1667        end do
1668 
1669      else
1670 
1671        write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
1672        do kk=0,nz-1
1673          do jj=0,ny-1
1674            do ii=0,nx-1
1675              ind=2*(1+ii+nx*(jj+ny*kk))-1
1676              write(denvtk,'(f16.8)') rhorfull(ind,1)
1677            end do
1678          end do
1679        end do
1680        write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double'
1681        do kk=0,nz-1
1682          do jj=0,ny-1
1683            do ii=0,nx-1
1684              ind=2*(1+ii+nx*(jj+ny*kk))
1685              write(denvtk,'(f16.8)') rhorfull(ind,1)
1686            end do
1687          end do
1688        end do
1689        write(denvtk,"(a,i18,a)") 'Re_mag 1 ',nfft_tot,' double'
1690        do kk=0,nz-1
1691          do jj=0,ny-1
1692            do ii=0,nx-1
1693              ind=2*(1+ii+nx*(jj+ny*kk))-1
1694              write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
1695            end do
1696          end do
1697        end do
1698        write(denvtk,"(a,i18,a)") 'Im_mag 1 ',nfft_tot,' double'
1699        do kk=0,nz-1
1700          do jj=0,ny-1
1701            do ii=0,nx-1
1702              ind=2*(1+ii+nx*(jj+ny*kk))
1703              write(denvtk,'(f16.8)') 2*rhorfull(ind,2)-rhorfull(ind,1)
1704            end do
1705          end do
1706        end do
1707 
1708      end if
1709 
1710    else  !here is the last option: nspden==4
1711 
1712      if(cplex==1) then
1713 
1714         !four component, write the density (scalar field) and magnetization density (vector field)
1715        write(denvtk,"(a,i18,a)") 'rho 1 ',nfft_tot,' double'
1716        do kk=0,nz-1
1717          do jj=0,ny-1
1718            do ii=0,nx-1
1719              ind=1+ii+nx*(jj+ny*kk)
1720              write(denvtk,'(f16.8)') rhorfull(ind,1)
1721            end do
1722          end do
1723        end do
1724        write(denvtk,"(a,i18,a)") 'mag 3 ',nfft_tot,' double'
1725        do kk=0,nz-1
1726          do jj=0,ny-1
1727            do ii=0,nx-1
1728              ind=1+ii+nx*(jj+ny*kk)
1729              write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
1730            end do
1731          end do
1732        end do
1733 
1734      else
1735 
1736        write(denvtk,"(a,i18,a)") 'Re_rho 1 ',nfft_tot,' double'
1737        do kk=0,nz-1
1738          do jj=0,ny-1
1739            do ii=0,nx-1
1740              ind=2*(1+ii+nx*(jj+ny*kk))-1
1741              write(denvtk,'(f16.8)') rhorfull(ind,1)
1742            end do
1743          end do
1744        end do
1745        write(denvtk,"(a,i18,a)") 'Im_rho 1 ',nfft_tot,' double'
1746        do kk=0,nz-1
1747          do jj=0,ny-1
1748            do ii=0,nx-1
1749              ind=2*(1+ii+nx*(jj+ny*kk))
1750              write(denvtk,'(f16.8)') rhorfull(ind,1)
1751            end do
1752          end do
1753        end do
1754        write(denvtk,"(a,i18,a)") 'Re_mag 3 ',nfft_tot,' double'
1755        do kk=0,nz-1
1756          do jj=0,ny-1
1757            do ii=0,nx-1
1758              ind=2*(1+ii+nx*(jj+ny*kk))-1
1759              write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
1760            end do
1761          end do
1762        end do
1763        write(denvtk,"(a,i18,a)") 'Im_mag 3 ',nfft_tot,' double'
1764        do kk=0,nz-1
1765          do jj=0,ny-1
1766            do ii=0,nx-1
1767              ind=2*(1+ii+nx*(jj+ny*kk))
1768              write(denvtk,'(3f16.8)') rhorfull(ind,2),rhorfull(ind,3),rhorfull(ind,4)
1769            end do
1770          end do
1771        end do
1772 
1773      end if
1774 
1775    end if ! nspden options condition
1776 
1777    close (denvtk)
1778 
1779     !clean up the gathered FFT mesh
1780    ABI_DEALLOCATE(rhorfull)
1781 
1782  end if
1783 
1784 end subroutine printmagvtk

m_dens/calcdensph [ Functions ]

[ Top ] [ m_dens ] [ Functions ]

NAME

 calcdensph

FUNCTION

 Compute and print integral of total density inside spheres around atoms.

INPUTS

  gmet(3,3)=reciprocal space metric tensor in bohr**-2
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell.
  nfft=(effective) number of FFT grid points (for this processor)
  ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
  nspden=number of spin-density components
  ntypat=number of atom types
  nunit=number of the unit for writing
  ratsph(ntypat)=radius of spheres around atoms
  rhor(nfft,nspden)=array for electron density in electrons/bohr**3.
   (total in first half and spin-up in second half if nspden=2)
   (total in first comp. and magnetization in comp. 2 to 4 if nspden=4)
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  typat(natom)=type of each atom
  ucvol=unit cell volume in bohr**3
  xred(3,natom)=reduced dimensionless atomic coordinates
  prtopt = if 1, the default printing is on, otherwise special printing options

OUTPUT

  intgden(nspden, natom)=intgrated density (magnetization...) for each atom in a sphere of radius ratsph. Optional arg
  dentot(nspden)=integrated density (magnetization...) over full u.c. vol, optional argument
  Rest is printing

PARENTS

      dfpt_scfcv,mag_constr,mag_constr_e,outscfcv

CHILDREN

      timab,wrtout,xmpi_sum

SOURCE

 864 subroutine calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,nunit,ratsph,rhor,rprimd,typat,ucvol,xred,&
 865 &    prtopt,cplex,intgden,dentot)
 866 
 867 
 868 !This section has been created automatically by the script Abilint (TD).
 869 !Do not modify the following lines by hand.
 870 #undef ABI_FUNC
 871 #define ABI_FUNC 'calcdensph'
 872 !End of the abilint section
 873 
 874  implicit none
 875 
 876 !Arguments ---------------------------------------------
 877 !scalars
 878  integer,intent(in)        :: natom,nfft,nspden,ntypat,nunit
 879  real(dp),intent(in)       :: ucvol
 880  type(MPI_type),intent(in) :: mpi_enreg
 881  integer ,intent(in)       :: prtopt
 882  integer, intent(in)       :: cplex
 883 !arrays
 884  integer,intent(in)  :: ngfft(18),typat(natom)
 885  real(dp),intent(in) :: gmet(3,3),ratsph(ntypat),rhor(cplex*nfft,nspden),rprimd(3,3)
 886  real(dp),intent(in) :: xred(3,natom)
 887 !integer,intent(out),optional   :: atgridpts(nfft)
 888  real(dp),intent(out),optional  :: intgden(nspden,natom)
 889  real(dp),intent(out),optional  :: dentot(nspden)
 890 !Local variables ------------------------------
 891 !scalars
 892  integer,parameter :: ishift=5
 893  integer :: i1,i2,i3,iatom,ierr,ifft_local,ix,iy,iz,izloc,n1,n1a,n1b,n2,ifft
 894  integer :: n2a,n2b,n3,n3a,n3b,nd3,nfftot
 895  integer :: ii
 896 !integer :: is,npts(natom)
 897  integer :: cmplex_den,jfft
 898  real(dp),parameter :: delta=0.99_dp
 899  real(dp) :: difx,dify,difz,r2,r2atsph,rr1,rr2,rr3,rx,ry,rz
 900  real(dp) :: fsm, ratsm, ratsm2
 901  real(dp) :: mag_coll   , mag_x, mag_y, mag_z ! EB
 902  real(dp) :: mag_coll_im, mag_x_im, mag_y_im, mag_z_im ! SPr
 903  real(dp) :: sum_mag, sum_mag_x,sum_mag_y,sum_mag_z,sum_rho_up,sum_rho_dn,sum_rho_tot ! EB
 904  real(dp) :: rho_tot, rho_tot_im
 905 ! real(dp) :: rho_up,rho_dn,rho_tot !EB
 906  logical   :: grid_found
 907  character(len=500) :: message
 908 !arrays
 909  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
 910  real(dp) :: intgden_(nspden,natom),tsec(2), my_xred(3, natom), xshift(3, natom)
 911 
 912 ! *************************************************************************
 913 
 914  n1=ngfft(1);n2=ngfft(2);n3=ngfft(3);nd3=n3/mpi_enreg%nproc_fft
 915  nfftot=n1*n2*n3
 916  intgden_=zero
 917 
 918  ! This routine is not able to handle xred positions that are "far" from the
 919  ! first unit cell so wrap xred into [0, 1[ interval here.
 920  call wrap2_zero_one(xred, my_xred, xshift)
 921 
 922  ratsm = zero
 923  if(present(intgden)) then
 924 
 925    ratsm = 0.05_dp ! default value for the smearing region radius - may become input variable later
 926 
 927  end if
 928 
 929 !Get the distrib associated with this fft_grid
 930  grid_found=.false.
 931 
 932  if(n2 == mpi_enreg%distribfft%n2_coarse ) then
 933 
 934    if(n3== size(mpi_enreg%distribfft%tab_fftdp3_distrib)) then
 935 
 936      fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3_distrib
 937      ffti3_local => mpi_enreg%distribfft%tab_fftdp3_local
 938      grid_found=.true.
 939 
 940    end if
 941 
 942  end if
 943 
 944  if(n2 == mpi_enreg%distribfft%n2_fine ) then
 945 
 946    if(n3 == size(mpi_enreg%distribfft%tab_fftdp3dg_distrib)) then
 947 
 948      fftn3_distrib => mpi_enreg%distribfft%tab_fftdp3dg_distrib
 949      ffti3_local => mpi_enreg%distribfft%tab_fftdp3dg_local
 950      grid_found = .true.
 951 
 952    end if
 953 
 954  end if
 955  if(.not.(grid_found)) then
 956 
 957    MSG_BUG("Unable to find an allocated distrib for this fft grid")
 958 
 959  end if
 960 
 961 !Loop over atoms
 962 !-------------------------------------------
 963  ii=0
 964 
 965  do iatom=1,natom
 966 
 967    !if (present(atgridpts)) then
 968    !  npts(iatom)=0           !SPr: initialize the number of grid points within atomic sphere around atom i
 969    !  ii=ii+1                 !SPr: initialize running index for constructing an array of grid point indexes
 970    !end if                    !     within atomic spheres
 971 
 972 !  Define a "box" around the atom
 973    r2atsph=1.0000001_dp*ratsph(typat(iatom))**2
 974    rr1=sqrt(r2atsph*gmet(1,1))
 975    rr2=sqrt(r2atsph*gmet(2,2))
 976    rr3=sqrt(r2atsph*gmet(3,3))
 977 
 978    n1a=int((my_xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1
 979    n1b=int((my_xred(1,iatom)+rr1+ishift)*n1      )-ishift*n1
 980    n2a=int((my_xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2
 981    n2b=int((my_xred(2,iatom)+rr2+ishift)*n2      )-ishift*n2
 982    n3a=int((my_xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3
 983    n3b=int((my_xred(3,iatom)+rr3+ishift)*n3      )-ishift*n3
 984 
 985    ratsm2 = (2*ratsph(typat(iatom))-ratsm)*ratsm
 986 
 987    do i3=n3a,n3b
 988      iz=mod(i3+ishift*n3,n3)
 989 
 990      if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
 991 
 992        izloc = ffti3_local(iz+1) - 1
 993        difz=dble(i3)/dble(n3)-my_xred(3,iatom)
 994        do i2=n2a,n2b
 995          iy=mod(i2+ishift*n2,n2)
 996          dify=dble(i2)/dble(n2)-my_xred(2,iatom)
 997          do i1=n1a,n1b
 998            ix=mod(i1+ishift*n1,n1)
 999 
1000            difx=dble(i1)/dble(n1)-my_xred(1,iatom)
1001            rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
1002            ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
1003            rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
1004            r2=rx**2+ry**2+rz**2
1005 
1006 
1007 !          Identify the fft indexes of the rectangular grid around the atom
1008            if(r2 > r2atsph) then
1009              cycle
1010            end if
1011 
1012            fsm = radsmear(r2, r2atsph, ratsm2)
1013 
1014            ifft_local=1+ix+n1*(iy+n2*izloc)
1015 
1016            !if(present(atgridpts)) then
1017            !  ii=ii+1
1018            !  atgridpts(ii)=ifft_local    !SPr: save grid point index (dbg: to check whether ifft_local is a valid "global" index )
1019            !end if
1020 
1021            if(nspden==1) then
1022 !            intgden_(1,iatom)= integral of total density
1023              intgden_(1,iatom)=intgden_(1,iatom)+fsm*rhor(ifft_local,1)
1024            elseif(nspden==2) then
1025 !            intgden_(1,iatom)= integral of up density
1026 !            intgden_(2,iatom)= integral of dn density
1027              intgden_(1,iatom)=intgden_(1,iatom)+fsm*rhor(ifft_local,2)
1028              intgden_(2,iatom)=intgden_(2,iatom)+fsm*rhor(ifft_local,1)-rhor(ifft_local,2)
1029            else
1030 !            intgden_(1,iatom)= integral of total density
1031 !            intgden_(2,iatom)= integral of magnetization, x-component
1032 !            intgden_(3,iatom)= integral of magnetization, y-component
1033 !            intgden_(4,iatom)= integral of magnetization, z-component
1034              intgden_(1,iatom)=intgden_(1,iatom)+fsm*rhor(ifft_local,1)
1035              intgden_(2,iatom)=intgden_(2,iatom)+fsm*rhor(ifft_local,2)
1036              intgden_(3,iatom)=intgden_(3,iatom)+fsm*rhor(ifft_local,3)
1037              intgden_(4,iatom)=intgden_(4,iatom)+fsm*rhor(ifft_local,4)
1038            end if
1039 
1040          end do
1041        end do
1042      end if
1043    end do
1044 
1045    intgden_(:,iatom)=intgden_(:,iatom)*ucvol/dble(nfftot)
1046 
1047    !if (present(atgridpts)) then
1048    !  npts(iatom)=ii-1
1049    !  do is=1,iatom-1,1
1050    !    npts(iatom)=npts(iatom)-npts(is)-1
1051    !  end do
1052    !  atgridpts(ii-npts(iatom))=npts(iatom)    !SPr: save number of grid points around atom i
1053    !end if
1054 
1055 
1056 !  End loop over atoms
1057 !  -------------------------------------------
1058 
1059  end do
1060 
1061 ! EB  - Compute magnetization of the whole cell
1062 ! SPr - in case of complex density array set cmplex_den to one
1063  cmplex_den=0
1064 
1065  if(cplex==2) then
1066 
1067    cmplex_den=1
1068 
1069  end if
1070 
1071  if(nspden==2) then
1072    mag_coll=zero
1073    mag_coll_im=zero
1074    rho_tot=zero
1075    rho_tot_im=zero
1076    do ifft=1,nfft
1077      jfft=(cmplex_den+1)*ifft
1078 !    rho_up=rho_up+rhor(ifft,2)
1079 !    rho_dn=rho_dn+(rhor(ifft,2)-rhor(ifft,1))
1080      rho_tot=rho_tot+rhor(jfft-cmplex_den,1)             ! real parts of density and magnetization
1081      mag_coll=mag_coll+2*rhor(jfft-cmplex_den,2)-rhor(jfft-cmplex_den,1)
1082 
1083      rho_tot_im=rho_tot_im+rhor(jfft,1)                  ! imaginary parts of density and magnetization
1084      mag_coll_im=mag_coll_im+2*rhor(jfft,2)-rhor(jfft,1) ! in case of real array both are equal to corresponding real parts
1085    end do
1086    mag_coll=mag_coll*ucvol/dble(nfftot)
1087    rho_tot =rho_tot *ucvol/dble(nfftot)
1088 
1089    mag_coll_im=mag_coll_im*ucvol/dble(nfftot)
1090    rho_tot_im =rho_tot_im *ucvol/dble(nfftot)
1091 !  rho_up=rho_up*ucvol/dble(nfftot)
1092 !  rho_dn=rho_dn*ucvol/dble(nfftot)
1093 !  rho_tot=rho_tot*ucvol/dble(nfftot)
1094  else if(nspden==4) then
1095    rho_tot=0
1096    rho_tot_im=0
1097    mag_x=0
1098    mag_y=0
1099    mag_z=0
1100    mag_x_im=0
1101    mag_y_im=0
1102    mag_z_im=0
1103    do ifft=1,nfft
1104      jfft=(cmplex_den+1)*ifft
1105      rho_tot=rho_tot+rhor(jfft-cmplex_den,1)
1106      mag_x=mag_x+rhor(jfft-cmplex_den,2)
1107      mag_y=mag_y+rhor(jfft-cmplex_den,3)
1108      mag_z=mag_z+rhor(jfft-cmplex_den,4)
1109      rho_tot_im=rho_tot_im+rhor(jfft,1)
1110      mag_x_im=mag_x_im+rhor(jfft,2)
1111      mag_y_im=mag_y_im+rhor(jfft,3)
1112      mag_z_im=mag_z_im+rhor(jfft,4)
1113    end do
1114    rho_tot=rho_tot*ucvol/dble(nfftot)
1115    mag_x=mag_x*ucvol/dble(nfftot)
1116    mag_y=mag_y*ucvol/dble(nfftot)
1117    mag_z=mag_z*ucvol/dble(nfftot)
1118    rho_tot_im=rho_tot_im*ucvol/dble(nfftot)
1119    mag_x_im=mag_x_im*ucvol/dble(nfftot)
1120    mag_y_im=mag_y_im*ucvol/dble(nfftot)
1121    mag_z_im=mag_z_im*ucvol/dble(nfftot)
1122  end if
1123 
1124 !MPI parallelization
1125  if(mpi_enreg%nproc_fft>1)then
1126    call timab(48,1,tsec)
1127    call xmpi_sum(intgden_,mpi_enreg%comm_fft,ierr)
1128    call xmpi_sum(rho_tot,mpi_enreg%comm_fft,ierr)  ! EB
1129    call xmpi_sum(mag_coll,mpi_enreg%comm_fft,ierr) ! EB
1130 
1131    call xmpi_sum(rho_tot_im,mpi_enreg%comm_fft,ierr)
1132    call xmpi_sum(mag_coll_im,mpi_enreg%comm_fft,ierr)
1133 !  call xmpi_sum(rho_up,mpi_enreg%comm_fft,ierr)  ! EB
1134 !  call xmpi_sum(rho_dn,mpi_enreg%comm_fft,ierr)  ! EB
1135    call xmpi_sum(mag_x,mpi_enreg%comm_fft,ierr)    ! EB
1136    call xmpi_sum(mag_y,mpi_enreg%comm_fft,ierr)    ! EB
1137    call xmpi_sum(mag_z,mpi_enreg%comm_fft,ierr)    ! EB
1138    call xmpi_sum(mag_x_im,mpi_enreg%comm_fft,ierr)    ! EB
1139    call xmpi_sum(mag_y_im,mpi_enreg%comm_fft,ierr)    ! EB
1140    call xmpi_sum(mag_z_im,mpi_enreg%comm_fft,ierr)    ! EB
1141    call timab(48,2,tsec)
1142  end if
1143 
1144 !Printing
1145  sum_mag=zero
1146  sum_mag_x=zero
1147  sum_mag_y=zero
1148  sum_mag_z=zero
1149  sum_rho_up=zero
1150  sum_rho_dn=zero
1151  sum_rho_tot=zero
1152 
1153  if(prtopt==1) then
1154 
1155    if(nspden==1) then
1156      write(message, '(4a)' ) &
1157 &     ' Integrated electronic density in atomic spheres:',ch10,&
1158 &     ' ------------------------------------------------'
1159      call wrtout(nunit,message,'COLL')
1160      write(message, '(a)' ) ' Atom  Sphere_radius  Integrated_density'
1161      call wrtout(nunit,message,'COLL')
1162      do iatom=1,natom
1163        write(message, '(i5,f15.5,f20.8)' ) iatom,ratsph(typat(iatom)),intgden_(1,iatom)
1164        call wrtout(nunit,message,'COLL')
1165      end do
1166    elseif(nspden==2) then
1167      write(message, '(4a)' ) &
1168 &     ' Integrated electronic and magnetization densities in atomic spheres:',ch10,&
1169 &     ' ---------------------------------------------------------------------'
1170      call wrtout(nunit,message,'COLL')
1171      write(message, '(3a)' ) ' Note: Diff(up-dn) is a rough ',&
1172 &     'approximation of local magnetic moment'
1173      call wrtout(nunit,message,'COLL')
1174      write(message, '(a)' ) ' Atom    Radius    up_density   dn_density  Total(up+dn)  Diff(up-dn)'
1175      call wrtout(nunit,message,'COLL')
1176      do iatom=1,natom
1177        write(message, '(i5,f10.5,2f13.6,a,f12.6,a,f12.6)' ) iatom,ratsph(typat(iatom)),intgden_(1,iatom),intgden_(2,iatom),&
1178 &       '  ',(intgden_(1,iatom)+intgden_(2,iatom)),' ',(intgden_(1,iatom)-intgden_(2,iatom))
1179        call wrtout(nunit,message,'COLL')
1180        ! Compute the sum of the magnetization
1181        sum_mag=sum_mag+intgden_(1,iatom)-intgden_(2,iatom)
1182        sum_rho_up=sum_rho_up+intgden_(1,iatom)
1183        sum_rho_dn=sum_rho_dn+intgden_(2,iatom)
1184        sum_rho_tot=sum_rho_tot+intgden_(1,iatom)+intgden_(2,iatom)
1185      end do
1186      write(message, '(a)') ' ---------------------------------------------------------------------'
1187      call wrtout(nunit,message,'COLL')
1188      write(message, '(a,2f13.6,a,f12.6,a,f12.6)') '  Sum:         ', sum_rho_up,sum_rho_dn,'  ',sum_rho_tot,' ',sum_mag
1189      call wrtout(nunit,message,'COLL')
1190      write(message, '(a,f14.6)') ' Total magnetization (from the atomic spheres):       ', sum_mag
1191      call wrtout(nunit,message,'COLL')
1192      write(message, '(a,f14.6)') ' Total magnetization (exact up - dn):                 ', mag_coll
1193      call wrtout(nunit,message,'COLL')
1194    ! EB for testing purpose print rho_up, rho_dn and rho_tot
1195 !    write(message, '(a,3f14.4,2i8)') ' rho_up, rho_dn, rho_tot, nfftot, nfft: ', rho_up,rho_dn,rho_tot,nfft,nfft
1196 !   call wrtout(nunit,message,'COLL')
1197 
1198    elseif(nspden==4) then
1199 
1200      write(message, '(4a)' ) &
1201 &     ' Integrated electronic and magnetization densities in atomic spheres:',ch10,&
1202 &     ' ---------------------------------------------------------------------'
1203      call wrtout(nunit,message,'COLL')
1204      write(message, '(3a)' ) ' Note:      this is a rough approximation of local magnetic moments'
1205      call wrtout(nunit,message,'COLL')
1206      write(message, '(a)' ) ' Atom   Radius      Total density     mag(x)      mag(y)      mag(z)  '
1207      call wrtout(nunit,message,'COLL')
1208      do iatom=1,natom
1209        write(message, '(i5,f10.5,f16.6,a,3f12.6)' ) iatom,ratsph(typat(iatom)),intgden_(1,iatom),'  ',(intgden_(ix,iatom),ix=2,4)
1210        call wrtout(nunit,message,'COLL')
1211        ! Compute the sum of the magnetization in x, y and z directions
1212        sum_mag_x=sum_mag_x+intgden_(2,iatom)
1213        sum_mag_y=sum_mag_y+intgden_(3,iatom)
1214        sum_mag_z=sum_mag_z+intgden_(4,iatom)
1215      end do
1216      write(message, '(a)') ' ---------------------------------------------------------------------'
1217      call wrtout(nunit,message,'COLL')
1218 !    write(message, '(a,f12.6,f12.6,f12.6)') ' Total magnetization :           ', sum_mag_x,sum_mag_y,sum_mag_z
1219      write(message, '(a,f12.6,f12.6,f12.6)') ' Total magnetization (spheres)   ', sum_mag_x,sum_mag_y,sum_mag_z
1220      call wrtout(nunit,message,'COLL')
1221      write(message, '(a,f12.6,f12.6,f12.6)') ' Total magnetization (exact)     ', mag_x,mag_y,mag_z
1222      call wrtout(nunit,message,'COLL')
1223 !    SPr for dfpt debug
1224 !    write(message, '(a,f12.6)') ' Total density (exact)           ', rho_tot
1225 !   call wrtout(nunit,message,'COLL')
1226    end if
1227 
1228  elseif(prtopt==-1) then
1229 
1230    write(message, '(2a)') ch10,' ------------------------------------------------------------------------'
1231    call wrtout(nunit,message,'COLL')
1232 
1233    if(nspden==1) then
1234      write(message, '(4a)' ) &
1235 &     ' Fermi level charge density n_f:',ch10,&
1236 &     ' ------------------------------------------------------------------------',ch10
1237    else
1238      write(message, '(4a)' ) &
1239 &     ' Fermi level charge density n_f and magnetization m_f:',ch10,&
1240 &     ' ------------------------------------------------------------------------',ch10
1241    end if
1242    call wrtout(nunit,message,'COLL')
1243 
1244    if(cmplex_den==0) then
1245      write(message, '(a,f13.8)') '     n_f   = ',rho_tot
1246    else
1247      write(message, '(a,f13.8,a,f13.8)') '  Re[n_f]= ', rho_tot,"   Im[n_f]= ",rho_tot_im
1248    end if
1249    call wrtout(nunit,message,'COLL')
1250    if(nspden==2) then
1251      if(cmplex_den==0) then
1252        write(message, '(a,f13.8)') '     m_f    = ', mag_coll
1253      else
1254        write(message, '(a,f13.8,a,f13.8)') '  Re[m_f]= ', mag_coll,"   Im[m_f]= ",mag_coll_im
1255      end if
1256      call wrtout(nunit,message,'COLL')
1257    elseif (nspden==4) then
1258      write(message, '(a,f13.8)') '     mx_f  = ',mag_x
1259      call wrtout(nunit,message,'COLL')
1260      write(message, '(a,f13.8)') '     my_f  = ',mag_y
1261      call wrtout(nunit,message,'COLL')
1262      write(message, '(a,f13.8)') '     mz_f  = ',mag_z
1263      call wrtout(nunit,message,'COLL')
1264    end if
1265 
1266    write(message, '(3a)') ch10,' ------------------------------------------------------------------------',ch10
1267    call wrtout(nunit,message,'COLL')
1268 
1269 
1270  else
1271  ! prtopt different from 1 and -1 (either 2,3 or 4)
1272 
1273    if(abs(rho_tot)<1.0d-10) then
1274      rho_tot=0
1275    end if
1276 
1277    write(message, '(2a)') ch10,' ------------------------------------------------------------------------'
1278    call wrtout(nunit,message,'COLL')
1279 
1280    if(nspden==1) then
1281      write(message, '(4a)' ) &
1282 &     ' Integral of the first order density n^(1):',ch10,&
1283 &     ' ------------------------------------------------------------------------',ch10
1284    else
1285      write(message, '(4a)' ) &
1286 &     ' Integrals of the first order density n^(1) and magnetization m^(1):',ch10,&
1287 &     ' ------------------------------------------------------------------------',ch10
1288    end if
1289    call wrtout(nunit,message,'COLL')
1290 
1291    if(cmplex_den==0) then
1292      write(message, '(a,e16.8)') '     n^(1)    = ', rho_tot
1293    else
1294      write(message, '(a,e16.8,a,e16.8)') '  Re[n^(1)] = ', rho_tot,"   Im[n^(1)] = ",rho_tot_im
1295    end if
1296    call wrtout(nunit,message,'COLL')
1297 
1298    if(nspden==2) then
1299 
1300      if(cmplex_den==0) then
1301        write(message, '(a,e16.8)') '     m^(1)    = ', mag_coll
1302      else
1303        write(message, '(a,e16.8,a,e16.8)') '  Re[m^(1)] = ', mag_coll,"   Im[m^(1)] = ",mag_coll_im
1304      end if
1305      call wrtout(nunit,message,'COLL')
1306 
1307    elseif (nspden==4) then
1308      if(cmplex_den==0) then
1309        write(message, '(a,e16.8)') '     mx^(1)   = ', mag_x
1310        call wrtout(nunit,message,'COLL')
1311        write(message, '(a,e16.8)') '     my^(1)   = ', mag_y
1312        call wrtout(nunit,message,'COLL')
1313        write(message, '(a,e16.8)') '     mz^(1)   = ', mag_z
1314        call wrtout(nunit,message,'COLL')
1315      else
1316        write(message, '(a,e16.8,a,e16.8)') '  Re[mx^(1)]= ',  mag_x, "   Im[mx^(1)]= ", mag_x_im
1317        call wrtout(nunit,message,'COLL')
1318        write(message, '(a,e16.8,a,e16.8)') '  Re[my^(1)]= ',  mag_y, "   Im[my^(1)]= ", mag_y_im
1319        call wrtout(nunit,message,'COLL')
1320        write(message, '(a,e16.8,a,e16.8)') '  Re[mz^(1)]= ',  mag_z, "   Im[mz^(1)]= ", mag_z_im
1321        call wrtout(nunit,message,'COLL')
1322      end if
1323    end if
1324 
1325    write(message, '(3a)') ch10,' ------------------------------------------------------------------------',ch10
1326    call wrtout(nunit,message,'COLL')
1327 
1328  end if
1329 
1330  if(present(intgden)) then
1331    intgden = intgden_
1332  end if
1333 
1334  if(present(dentot)) then
1335    if(nspden==2) then
1336      dentot(1)=rho_tot
1337      dentot(2)=mag_coll
1338    elseif(nspden==4) then
1339      dentot(1)=rho_tot
1340      dentot(2)=mag_x
1341      dentot(3)=mag_y
1342      dentot(4)=mag_z
1343    end if
1344  end if
1345 
1346 end subroutine calcdensph

m_dens/dens_hirsh [ Functions ]

[ Top ] [ m_dens ] [ Functions ]

NAME

 dens_hirsh

FUNCTION

 Compute the Hirshfeld charges

INPUTS

  mpoint=Maximum number of points in radial meshes.
  radii(mpoint, ntypat)=Radial meshes for each type
  aeden(mpoint, nytpat)=All-electron densities.
  npoint(ntypat)=The number of the last point with significant density is stored in npoint(itypat)
  minimal_den=Tolerance on the minum value of the density
  grid_den(nrx,nry,nrz)= density on the grid
  natom = number of atoms in the unit cell
  nrx,nry,nrz= number of points in the grid for the three directions
  ntypat=number of types of atoms in unit cell.
  rprimd(3,3)=dimensional primitive translations in real space (bohr)
  typat(natom)=type of each atom
  xcart(3,natom) = different positions of the atoms in the unit cell
  zion=(ntypat)gives the ionic charge for each type of atom
  prtcharge=1 to write the Hirshfeld charge decomposition

OUTPUT

  hcharge(natom), hden(natom), hweight(natom)= Hirshfeld charges, densities, weights.

PARENTS

      m_cut3d

CHILDREN

      metric,spline,xcart2xred

SOURCE

 89 subroutine dens_hirsh(mpoint,radii,aeden,npoint,minimal_den,grid_den, &
 90   natom,nrx,nry,nrz,ntypat,rprimd,xcart,typat,zion,prtcharge,hcharge,hden,hweight)
 91 
 92 
 93 !This section has been created automatically by the script Abilint (TD).
 94 !Do not modify the following lines by hand.
 95 #undef ABI_FUNC
 96 #define ABI_FUNC 'dens_hirsh'
 97 !End of the abilint section
 98 
 99  implicit none
100 
101 !Arguments ------------------------------------
102 !scalars
103  integer,intent(in) :: natom,nrx,nry,nrz,ntypat,prtcharge,mpoint
104  real(dp),intent(in) :: minimal_den
105 !arrays
106  integer,intent(in) :: typat(natom),npoint(ntypat)
107  real(dp),intent(in) :: grid_den(nrx,nry,nrz),rprimd(3,3),zion(ntypat)
108  real(dp),intent(in) :: xcart(3,natom)
109  real(dp),intent(in) :: radii(mpoint,ntypat),aeden(mpoint,ntypat)
110  real(dp),intent(out) :: hcharge(natom),hden(natom),hweight(natom)
111 
112 !Local variables -------------------------
113 !scalars
114  integer :: i1,i2,i3,iatom,icell,igrid,ii,inmax,inmin,istep,itypat
115  integer :: k1,k2,k3,mcells,nfftot,ngoodpoints,npt
116  real(dp) :: aa,bb,coeff1,coeff2,coeff3,den,factor,h_inv,hh,maxrad
117  real(dp) :: rr,rr2,total_charge,total_weight,total_zion,ucvol
118  real(dp) :: yp1,ypn
119  !character(len=500) :: msg
120 !arrays
121  integer :: highest(3),lowest(3)
122  integer,allocatable :: ncells(:)
123  real(dp) :: coordat(3),coord23_1,coord23_2,coord23_3,diff1,diff2,diff3,gmet(3,3),gprimd(3,3),rmet(3,3)
124  real(dp) :: vperp(3),width(3)
125  real(dp),allocatable :: coord1(:,:),local_den(:,:,:,:)
126  real(dp),allocatable :: step(:,:),sum_den(:,:,:),work(:)
127  real(dp),allocatable :: xcartcells(:,:,:),xred(:,:),yder2(:)
128 
129 ! *********************************************************************
130 
131 !1. Read the 1D all-electron atomic files
132 !Store the radii in radii(:,itypat), and the all-electron
133 !densities in aeden(:,itypat). The number of the last
134 !point with significant density is stored in npoint(itypat)
135 
136 !2. Compute the list of atoms that are sufficiently close to the cell
137 
138  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
139  nfftot=nrx*nry*nrz
140 !DEBUG
141 !do k3=1,nrz
142 !do k2=1,nry
143 !do k1=1,nrx
144 !total_charge=total_charge+grid_den(k1,k2,k3)
145 !end do
146 !end do
147 !end do
148 !write(std_out,*)' total_charge=',total_charge*ucvol/dble(nfftot)
149 !ENDDEBUG
150 
151  ABI_ALLOCATE(xred,(3,natom))
152  call xcart2xred(natom,rprimd,xcart,xred)
153 
154 !Compute the widths of the cell
155 !First width : perpendicular vector length
156  vperp(:)=rprimd(:,1)-rprimd(:,2)*rmet(1,2)/rmet(2,2) -rprimd(:,3)*rmet(1,3)/rmet(3,3)
157  width(1)=sqrt(dot_product(vperp,vperp))
158 !Second width
159  vperp(:)=rprimd(:,2)-rprimd(:,1)*rmet(2,1)/rmet(1,1) -rprimd(:,3)*rmet(2,3)/rmet(3,3)
160  width(2)=sqrt(dot_product(vperp,vperp))
161 !Third width
162  vperp(:)=rprimd(:,3)-rprimd(:,1)*rmet(3,1)/rmet(1,1) -rprimd(:,2)*rmet(3,2)/rmet(2,2)
163  width(3)=sqrt(dot_product(vperp,vperp))
164 
165 !Compute the number of cells that will make up the supercell
166  ABI_ALLOCATE(ncells,(natom))
167  mcells=0
168  do iatom=1,natom
169    itypat=typat(iatom)
170    maxrad=radii(npoint(itypat),itypat)
171 !  Compute the lower and higher indices of the supercell
172 !  for this atom
173    do ii=1,3
174      lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii))
175      highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1)
176 !    Next coding, still incorrect
177 !    lowest(ii)=floor(xred(ii,iatom)-maxrad/width(ii))-1
178 !    highest(ii)=ceiling(xred(ii,iatom)+maxrad/width(ii))+1
179 !    Old coding incorrect
180 !    lowest(ii)=ceiling(-xred(ii,iatom)-maxrad/width(ii))
181 !    highest(ii)=floor(-xred(ii,iatom)+maxrad/width(ii)+1)
182    end do
183    ncells(iatom)=(highest(1)-lowest(1)+1)* &
184 &   (highest(2)-lowest(2)+1)* &
185 &   (highest(3)-lowest(3)+1)
186 !  DEBUG
187 !  write(std_out,*)' maxrad=',maxrad
188 !  write(std_out,*)' lowest(:)=',lowest(:)
189 !  write(std_out,*)' highest(:)=',highest(:)
190 !  write(std_out,*)' ncells(iatom)=',ncells(iatom)
191 !  ENDDEBUG
192  end do
193  mcells=maxval(ncells(:))
194 
195 !Compute, for each atom, the set of image atoms in the whole supercell
196  ABI_ALLOCATE(xcartcells,(3,mcells,natom))
197  do iatom=1,natom
198    itypat=typat(iatom)
199    maxrad=radii(npoint(itypat),itypat)
200 !  Compute the lower and higher indices of the supercell
201 !  for this atom
202 
203    do ii=1,3
204      lowest(ii)=floor(-xred(ii,iatom)-maxrad/width(ii))
205      highest(ii)=ceiling(-xred(ii,iatom)+maxrad/width(ii)+1)
206    end do
207    icell=0
208    do i1=lowest(1),highest(1)
209      do i2=lowest(2),highest(2)
210        do i3=lowest(3),highest(3)
211          icell=icell+1
212          xcartcells(:,icell,iatom)=xcart(:,iatom)+i1*rprimd(:,1)+i2*rprimd(:,2)+i3*rprimd(:,3)
213        end do
214      end do
215    end do
216  end do
217 
218 !Compute, for each atom, the all-electron pro-atom density
219 !at each point in the primitive cell
220  ABI_ALLOCATE(local_den,(nrx,nry,nrz,natom))
221  ABI_ALLOCATE(step,(2,mpoint))
222  ABI_ALLOCATE(work,(mpoint))
223  ABI_ALLOCATE(yder2,(mpoint))
224  ABI_ALLOCATE(coord1,(3,nrx))
225  coeff1=one/nrx
226  coeff2=one/nry
227  coeff3=one/nrz
228 
229  do iatom=1,natom
230    itypat=typat(iatom)
231    npt=npoint(itypat)
232    maxrad=radii(npt,itypat)
233 !   write(std_out,*)
234 !   write(std_out,'(a,i4)' )' hirsh : accumulating density for atom ',iatom
235 !  write(std_out,*)' ncells(iatom)=',ncells(iatom)
236    do istep=1,npt-1
237      step(1,istep)=radii(istep+1,itypat) - radii(istep,itypat)
238      step(2,istep)=one/step(1,istep)
239    end do
240 !  Approximate first derivative for small radii
241    yp1=(aeden(2,itypat)-aeden(1,itypat))/(radii(2,itypat)-radii(1,itypat))
242    ypn=zero
243    call spline(radii(1:npt,itypat),aeden(1:npt,itypat),npt,yp1,ypn,yder2)
244 
245    local_den(:,:,:,iatom)=zero
246 
247 !  Big loop on the cells
248    do icell=1,ncells(iatom)
249 !    write(std_out,*)' icell=',icell
250      coordat(:)=xcartcells(:,icell,iatom)
251 
252 !    Big loop on the grid points
253      do k1 = 1,nrx
254        coord1(:,k1)=rprimd(:,1)*(k1-1)*coeff1
255      end do
256      do k3 = 1, nrz
257        do k2 = 1, nry
258          coord23_1=rprimd(1,2)*(k2-1)*coeff2+rprimd(1,3)*(k3-1)*coeff3-coordat(1)
259          coord23_2=rprimd(2,2)*(k2-1)*coeff2+rprimd(2,3)*(k3-1)*coeff3-coordat(2)
260          coord23_3=rprimd(3,2)*(k2-1)*coeff2+rprimd(3,3)*(k3-1)*coeff3-coordat(3)
261          do k1 = 1, nrx
262            diff1=coord1(1,k1)+coord23_1
263            diff2=coord1(2,k1)+coord23_2
264            diff3=coord1(3,k1)+coord23_3
265            rr2=diff1**2+diff2**2+diff3**2
266            if(rr2<maxrad**2)then
267 
268              rr=sqrt(rr2)
269 !            Find the index of the radius by bissection
270              if (rr < radii(1,itypat)) then
271 !              Linear extrapolation
272                den=aeden(1,itypat)+(rr-radii(1,itypat))/(radii(2,itypat)-radii(1,itypat))&
273 &               *(aeden(2,itypat)-aeden(1,itypat))
274              else
275 !              Use the spline interpolation
276 !              Find the index of the radius by bissection
277                inmin=1
278                inmax=npt
279                igrid=1
280                do
281                  if(inmax-inmin==1)exit
282                  igrid=(inmin+inmax)/2
283                  if(rr>=radii(igrid,itypat))then
284                    inmin=igrid
285                  else
286                    inmax=igrid
287                  end if
288                end do
289                igrid=inmin
290 !              write(std_out,*)' igrid',igrid
291 
292                hh=step(1,igrid)
293                h_inv=step(2,igrid)
294                aa= (radii(igrid+1,itypat)-rr)*h_inv
295                bb= (rr-radii(igrid,itypat))*h_inv
296                den = aa*aeden(igrid,itypat) + bb*aeden(igrid+1,itypat)  &
297 &               +( (aa*aa*aa-aa)*yder2(igrid)         &
298 &               +(bb*bb*bb-bb)*yder2(igrid+1) ) *hh*hh*sixth
299              end if ! Select small radius or spline
300 
301              local_den(k1,k2,k3,iatom)=local_den(k1,k2,k3,iatom)+den
302            end if ! dist2<maxrad
303 
304          end do ! k1
305        end do ! k2
306      end do ! k3
307 
308    end do ! icell
309  end do ! iatom
310 
311 !Compute, the total all-electron density at each point in the primitive cell
312  ABI_ALLOCATE(sum_den,(nrx,nry,nrz))
313  sum_den(:,:,:)=zero
314  do iatom=1,natom
315    sum_den(:,:,:)=sum_den(:,:,:)+local_den(:,:,:,iatom)
316  end do
317 
318 !DEBUG
319 !do k3=1,nrz
320 !do k2=1,nry
321 !do k1=1,nrx
322 !write(std_out,'(3i4,3es16.6)' )k1,k2,k3,local_den(k1,k2,k3,1:2),sum_den(k1,k2,k3)
323 !end do
324 !end do
325 !end do
326 !write(std_out,*)' hirsh : before accumulate the integral of the density'
327 !ENDDEBUG
328 
329 
330 !Accumulate the integral of the density, to get Hirshfeld charges
331 !There is a minus sign because the electron has a negative charge
332  ngoodpoints = 0
333  hcharge(:)=zero
334  hweight(:)=zero
335  do k3=1,nrz
336    do k2=1,nry
337      do k1=1,nrx
338 !      Use minimal_den in order to avoid divide by zero
339        if (abs(sum_den(k1,k2,k3)) > minimal_den) then
340          ngoodpoints = ngoodpoints+1
341          factor=grid_den(k1,k2,k3)/(sum_den(k1,k2,k3)+minimal_den)
342          do iatom=1,natom
343            hden(iatom)=hden(iatom)+local_den(k1,k2,k3,iatom)
344            hcharge(iatom)=hcharge(iatom)-local_den(k1,k2,k3,iatom)*factor
345            hweight(iatom)=hweight(iatom)+local_den(k1,k2,k3,iatom)/(sum_den(k1,k2,k3)+minimal_den)
346          end do
347        end if
348      end do
349    end do
350  end do
351 
352 !DEBUG
353 !do iatom=1,natom
354 !write(std_out,'(i9,3es17.6)' )iatom,hden(iatom),hcharge(iatom),hweight(iatom)
355 !end do
356 !ENDDEBUG
357 
358  hcharge(:)=hcharge(:)*ucvol/dble(nfftot)
359  hweight(:)=hweight(:)/dble(nfftot)
360 
361 !Check on the total charge
362  total_zion=sum(zion(typat(1:natom)))
363  total_charge=sum(hcharge(1:natom))
364  total_weight=sum(hweight(1:natom))
365 
366 !DEBUG
367 !write(std_out,*)' ngoodpoints = ', ngoodpoints, ' out of ', nfftot
368 !write(std_out,*)' total_weight=',total_weight
369 !write(std_out,*)' total_weight=',total_weight
370 !ENDDEBUG
371 
372 !Output
373  if (prtcharge == 1) then
374      write(std_out,*)
375      write(std_out,*)'    Hirshfeld analysis'
376      write(std_out,*)'    Atom       Zion       Electron  Charge       Net charge '
377      write(std_out,*)
378      do iatom=1,natom
379        write(std_out,'(i9,3es17.6)' )&
380 &       iatom,zion(typat(iatom)),hcharge(iatom),hcharge(iatom)+zion(typat(iatom))
381      end do
382      write(std_out,*)
383      write(std_out,'(a,3es17.6)')'    Total',total_zion,total_charge,total_charge+total_zion
384      write(std_out,*)
385  end if
386 
387  ABI_DEALLOCATE(coord1)
388  ABI_DEALLOCATE(local_den)
389  ABI_DEALLOCATE(ncells)
390  ABI_DEALLOCATE(step)
391  ABI_DEALLOCATE(sum_den)
392  ABI_DEALLOCATE(work)
393  ABI_DEALLOCATE(xcartcells)
394  ABI_DEALLOCATE(xred)
395  ABI_DEALLOCATE(yder2)
396 
397 end subroutine dens_hirsh

m_dens/mag_constr [ Functions ]

[ Top ] [ m_dens ] [ Functions ]

NAME

 mag_constr

FUNCTION

 This routine is called to compute the potential corresponding to constrained magnetic moments.

INPUTS

  natom=number of atoms
  spinat=fixed magnetic moments vectors
  nspden = number of spin densities (1 2 or 4)
  magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size
  magcon_lambda=the size of the penalty terms
  rprimd=lattice vectors (dimensionful)
  mpi_enreg=mpi structure with communicator info
  nfft=number of points in standard fft grid
  ngfft=FFT grid dimensions
  ntypat=number of types of atoms
  ratsph=radii for muffin tin spheres of each atom
  rhor=density in real space
  typat=types of atoms
  xred=reduced atomic positions

OUTPUT

  Vmagconstr=the constraining potential

PARENTS

      energy,rhotov,setvtr

CHILDREN

      calcdensph,metric,ptabs_fourdp,timab,xmpi_sum

NOTES

  based on html notes for the VASP implementation at
  http://cms.mpi.univie.ac.at/vasp/vasp/Constraining_direction_magnetic_moments.html

SOURCE

438 subroutine mag_constr(natom,spinat,nspden,magconon,magcon_lambda,rprimd, &
439                       mpi_enreg,nfft,ngfft,ntypat,ratsph,rhor, &
440                       typat,Vmagconstr,xred)
441 
442 
443 !This section has been created automatically by the script Abilint (TD).
444 !Do not modify the following lines by hand.
445 #undef ABI_FUNC
446 #define ABI_FUNC 'mag_constr'
447 !End of the abilint section
448 
449  implicit none
450 
451 !Arguments ------------------------------------
452 !scalars
453  integer,intent(in) :: natom,magconon,nfft,nspden
454  integer,intent(in) :: ntypat
455  real(dp),intent(in) :: magcon_lambda
456  real(dp),intent(out) :: Vmagconstr(nfft,nspden)
457  type(MPI_type),intent(in) :: mpi_enreg
458 !arrays
459  integer,intent(in)  :: typat(natom)
460  integer,intent(in)  :: ngfft(18)
461  real(dp),intent(in) :: ratsph(ntypat)
462  real(dp),intent(in) :: rhor(nfft,nspden)
463  real(dp),intent(in) :: rprimd(3,3)
464  real(dp),intent(in) :: spinat(3,natom)
465  real(dp),intent(in) :: xred(3,natom)
466 
467 !Local variables-------------------------------
468 !scalars
469  integer,parameter :: ishift=5
470  integer :: iatom, ierr
471  integer :: cplex1=1
472  integer :: n1a, n1b, n3a, n3b, n2a, n2b
473  integer :: n1, n2, n3
474  integer :: ifft_local
475  integer ::  i1,i2,i3,ix,iy,iz,izloc,nd3
476  real(dp) :: arg,intgden_proj,r2atsph,rr1,rr2,rr3,fsm,ratsm,ratsm2,difx,dify,difz,r2,rx,ry,rz
477  real(dp) :: norm
478  real(dp),parameter :: delta=0.99_dp
479 !arrays
480  real(dp) :: intgden(nspden,natom)
481  real(dp) :: spinat_norm(3,natom)
482  !real(dp),allocatable :: cmm_x(:,:),cmm_y(:,:),cmm_z(:,:)
483  real(dp):: cmm_x,cmm_y,cmm_z
484  real(dp) :: gprimd(3,3),rmet(3,3),gmet(3,3),ucvol
485  real(dp) :: tsec(2)
486  integer, ABI_CONTIGUOUS pointer :: fftn2_distrib(:),ffti2_local(:)
487  integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
488 
489 ! ***********************************************************************************************
490 
491 !We need the metric because it is needed in calcdensph.F90
492  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
493 
494 !We need the integrated magnetic moments and the smoothing function
495  call calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,std_out,ratsph,rhor,rprimd,typat,ucvol,xred,1,cplex1,intgden=intgden)
496 
497  n1 = ngfft(1)
498  n2 = ngfft(2)
499  n3 = ngfft(3)
500  nd3=n3/mpi_enreg%nproc_fft
501 
502  ratsm = 0.05_dp ! default value for the smearing region radius - may become input variable later
503  Vmagconstr = zero
504 
505 !Get the distrib associated with this fft_grid
506  call ptabs_fourdp(mpi_enreg,n2,n3,fftn2_distrib,ffti2_local,fftn3_distrib,ffti3_local)
507 
508 !Loop over atoms
509 !-------------------------------------------
510  do iatom=1,natom
511 
512    norm = sqrt(sum(spinat(:,iatom)**2))
513    spinat_norm(:,iatom) = zero
514    if (norm > tol10) then
515      spinat_norm(:,iatom) = spinat(:,iatom) / norm
516    else if (magconon == 1) then
517 !    if spinat = 0 and we are imposing the direction only, skip this atom
518      cycle
519    end if
520 
521 !  Define a "box" around the atom
522    r2atsph=1.0000001_dp*ratsph(typat(iatom))**2
523    rr1=sqrt(r2atsph*gmet(1,1))
524    rr2=sqrt(r2atsph*gmet(2,2))
525    rr3=sqrt(r2atsph*gmet(3,3))
526 
527    n1a=int((xred(1,iatom)-rr1+ishift)*n1+delta)-ishift*n1
528    n1b=int((xred(1,iatom)+rr1+ishift)*n1      )-ishift*n1
529    n2a=int((xred(2,iatom)-rr2+ishift)*n2+delta)-ishift*n2
530    n2b=int((xred(2,iatom)+rr2+ishift)*n2      )-ishift*n2
531    n3a=int((xred(3,iatom)-rr3+ishift)*n3+delta)-ishift*n3
532    n3b=int((xred(3,iatom)+rr3+ishift)*n3      )-ishift*n3
533 
534    ratsm2 = -(ratsm**2 - 2*ratsph(typat(iatom))*ratsm)
535 
536    do i3=n3a,n3b
537      iz=mod(i3+ishift*n3,n3)
538      if(fftn3_distrib(iz+1)==mpi_enreg%me_fft) then
539        izloc = ffti3_local(iz+1) - 1
540        difz=dble(i3)/dble(n3)-xred(3,iatom)
541        do i2=n2a,n2b
542          iy=mod(i2+ishift*n2,n2)
543          dify=dble(i2)/dble(n2)-xred(2,iatom)
544          do i1=n1a,n1b
545            ix=mod(i1+ishift*n1,n1)
546            difx=dble(i1)/dble(n1)-xred(1,iatom)
547            rx=difx*rprimd(1,1)+dify*rprimd(1,2)+difz*rprimd(1,3)
548            ry=difx*rprimd(2,1)+dify*rprimd(2,2)+difz*rprimd(2,3)
549            rz=difx*rprimd(3,1)+dify*rprimd(3,2)+difz*rprimd(3,3)
550            r2=rx**2+ry**2+rz**2
551 
552 
553 !          Identify the fft indexes of the rectangular grid around the atom
554            if (r2 > r2atsph) cycle
555 
556            fsm = radsmear(r2, r2atsph, ratsm2)
557 
558            ifft_local=1+ix+n1*(iy+n2*izloc)
559 
560 !          Calculate the x- and y-components of the square bracket term
561            cmm_x = zero
562            cmm_y = zero
563            cmm_z = zero
564            intgden_proj = zero
565            if (nspden == 4) then
566              if (magconon==1) then
567 !              Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector
568 !              This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector
569                intgden_proj=spinat_norm(1,iatom)*intgden(2,iatom)+ &
570 &               spinat_norm(2,iatom)*intgden(3,iatom)+ &
571 &               spinat_norm(3,iatom)*intgden(4,iatom)
572 
573                cmm_x=intgden(2,iatom)
574                cmm_x=cmm_x-spinat_norm(1,iatom)*intgden_proj
575 
576                cmm_y=intgden(3,iatom)
577                cmm_y=cmm_y-spinat_norm(2,iatom)*intgden_proj
578 
579              else if (magconon==2 .and. nspden == 4) then
580                cmm_x=intgden(2,iatom)-spinat(1,iatom)
581                cmm_y=intgden(3,iatom)-spinat(2,iatom)
582              end if
583 
584 !            do loop on spirals - should be added later on. arg = scalar product k.R (with 2pi factor?)
585              arg=0
586 
587 !            Calculate the constraining potential for x- and y- components of the mag. mom. vector
588 !            Eric Bousquet has derived the relationship between spin components and potential spin matrix elements:
589 !            1 = up up     = +z
590 !            2 = down down = -z
591 !            3 = up down   = +x
592 !            4 = down up   = -y
593              Vmagconstr(ifft_local,3)=Vmagconstr(ifft_local,3) + &
594 &             2*magcon_lambda*fsm*cmm_x
595 !            & 2*magcon_lambda*fsm*(cmm_x*cos(arg)+cmm_y*sin(arg))
596              Vmagconstr(ifft_local,4)=Vmagconstr(ifft_local,4) - &
597 &             2*magcon_lambda*fsm*cmm_y
598 !            & 2*magcon_lambda*fsm*(cmm_y*cos(arg)+cmm_x*sin(arg))
599 !            end do loop on spirals
600            end if ! nspden 4
601 
602 !          Calculate the z-component of the square bracket term
603            if (magconon==1) then
604              if (nspden == 4) then
605                ! m_z - spinat_z * <m | spinat>
606                cmm_z = intgden(4,iatom) - spinat_norm(3,iatom)*intgden_proj
607              else if (nspden == 2) then
608                ! this will be just a sign +/- : are we in the same direction as spinat_z?
609                !    need something more continuous??? To make sure the gradient pushes the state towards FM/AFM?
610                cmm_z = -sign(one, (intgden(1,iatom)-intgden(2,iatom))*spinat_norm(3,iatom))
611              end if
612            else if (magconon==2) then
613              if (nspden == 4) then
614                cmm_z=intgden(4,iatom)-spinat(3,iatom)
615              else if (nspden == 2) then
616                ! this is up spins - down spins - requested moment ~ 0
617                ! EB: note that intgden comes from calcdensph, which, in nspden=2 case, returns
618                ! intgden(1)=rho_up=n+m
619                ! intgden(2)=rho_dn=n-m
620                ! Then, is the following line be
621                ! cmm_z=half*(intgden(1,iatom)-intgden(2,iatom)) - spinat(3,iatom)
622                ! ??
623                cmm_z=intgden(1,iatom)-intgden(2,iatom) - spinat(3,iatom)
624              end if
625            end if
626 
627 !          Calculate the constraining potential for z-component of the mag. mom. vector
628            Vmagconstr(ifft_local,1)=Vmagconstr(ifft_local,1) + 2*magcon_lambda*fsm*cmm_z
629            Vmagconstr(ifft_local,2)=Vmagconstr(ifft_local,2) - 2*magcon_lambda*fsm*cmm_z
630 !          end do loop on spirals
631 
632          end do  ! i1
633        end do  ! i2
634      end if  ! if this is my fft slice
635    end do ! i3
636 
637 !  end loop over atoms
638  end do
639 
640 !MPI parallelization
641 !TODO: test if xmpi_sum does the correct stuff for a slice of Vmagconstr
642  if(mpi_enreg%nproc_fft>1)then
643    call timab(48,1,tsec)
644    call xmpi_sum(Vmagconstr,mpi_enreg%comm_fft,ierr)
645    call timab(48,2,tsec)
646  end if
647 
648 ! write (201,*) '# potential 1'
649 ! write (201,*) Vmagconstr(:,1)
650 
651 ! write (202,*) '# potential 2'
652 ! write (202,*) Vmagconstr(:,2)
653 
654 ! if (nspden > 2) then
655 !   write (203,*) '# potential 3'
656 !   write (203,*) Vmagconstr(:,3)
657 
658 !   write (204,*) '# potential 4'
659 !   write (204,*) Vmagconstr(:,4)
660 ! end if
661 
662 end subroutine mag_constr

m_dens/mag_constr_e [ Functions ]

[ Top ] [ m_dens ] [ Functions ]

NAME

 mag_constr_e

FUNCTION

 Compute the energy corresponding to constrained magnetic moments.

INPUTS

  magconon=constraining option (on/off); 1=fix only the direction, 2=fix the direction and size
  spinat=fixed magnetic moments vectors
  magcon_lambda=the size of the penalty terms

OUTPUT

  Epen=penalty contribution to the total energy corresponding to the constrained potential
  Econstr=???
  Eexp=???

PARENTS

      outscfcv

CHILDREN

      calcdensph,metric,wrtout

SOURCE

691 subroutine mag_constr_e(magconon,magcon_lambda,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,ratsph,rhor,rprimd,spinat,typat,xred)
692 
693 
694 !This section has been created automatically by the script Abilint (TD).
695 !Do not modify the following lines by hand.
696 #undef ABI_FUNC
697 #define ABI_FUNC 'mag_constr_e'
698 !End of the abilint section
699 
700  implicit none
701 
702 !Arguments ------------------------------------
703 !scalars
704  integer,intent(in) :: natom,magconon,nspden,nfft,ntypat
705  real(dp),intent(in) :: magcon_lambda
706 !arrays
707  integer, intent(in) :: ngfft(18),typat(natom)
708  real(dp),intent(in) :: spinat(3,natom), rprimd(3,3)
709  real(dp),intent(in) :: ratsph(ntypat),rhor(nfft,nspden),xred(3,natom)
710  type(MPI_type),intent(in) :: mpi_enreg
711 
712 !Local variables-------------------------------
713 !scalars
714  integer :: iatom,ii
715  integer :: cplex1=1    ! dummy argument for calcdensphere
716  real(dp) :: intgden_proj, Epen,Econstr,lVp, norm
717 !arrays
718  real(dp) :: intmm(3), mag_1atom(3)
719  real(dp), allocatable :: intgden(:,:)
720  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),ucvol
721  real(dp) :: spinat_norm(3,natom)
722  character(len=500) :: message
723 
724 ! *********************************************************************
725 
726 !We need the metric because it is needed in calcdensph.F90
727  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
728 
729  ABI_ALLOCATE (intgden, (nspden,natom))
730 
731 !We need the integrated magnetic moments
732  cplex1=1
733  call calcdensph(gmet,mpi_enreg,natom,nfft,ngfft,nspden,ntypat,std_out,ratsph,rhor,rprimd,typat,ucvol,xred,&
734 & 1,cplex1,intgden=intgden)
735 
736  Epen=0
737  Econstr=0
738  lVp=0
739 
740 !Loop over atoms
741 !-------------------------------------------
742  do iatom=1,natom
743 
744    norm = sqrt(sum(spinat(:,iatom)**2))
745    spinat_norm(:,iatom) = zero
746    if (norm > tol10) then
747      spinat_norm(:,iatom) = spinat(:,iatom) / norm
748    else if (magconon == 1) then
749 !    if spinat = 0 and we are imposing the direction only, skip this atom
750      cycle
751    end if
752 !  Calculate the scalar product of the fixed mag. mom. vector and calculated mag. mom. vector
753 !  This is actually the size of the projection of the calc. mag. mom. vector on the fixed mag. mom. vector
754 
755 ! for the collinear spin case, set up a fictitious 3D vector along z
756    if (nspden == 4) then
757      mag_1atom(1:3) = intgden(2:4,iatom)
758    else if (nspden == 2) then
759      mag_1atom = zero
760      mag_1atom(3) = intgden(1,iatom)-intgden(2,iatom)
761    end if
762 
763    intgden_proj = zero
764    intmm = zero
765 !  Calculate the square bracket term
766    if (magconon==1) then
767      intgden_proj=spinat_norm(1,iatom)*mag_1atom(1)+ &
768 &     spinat_norm(2,iatom)*mag_1atom(2)+ &
769 &     spinat_norm(3,iatom)*mag_1atom(3)
770 
771      do ii=1,3
772        intmm(ii)=mag_1atom(ii)-spinat_norm(ii,iatom)*intgden_proj
773      end do
774 
775 !    Calculate the energy Epen corresponding to the constraining potential
776 !    Econstr and lVp do not have a clear meaning (yet)
777      Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3))
778      Econstr=Econstr-magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3))
779      lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3))
780 
781    else if (magconon==2) then
782      do ii=1,3
783        intmm(ii)=mag_1atom(ii)-spinat(ii,iatom)
784      end do
785 
786 !    Calculate the energy Epen corresponding to the constraining potential
787 !    Epen = -Econstr - lVp
788 !    Econstr = -M**2 + spinat**2
789 !    lVp = +2 M \cdot spinat
790      Epen=Epen+magcon_lambda*(intmm(1)*intmm(1)+intmm(2)*intmm(2)+intmm(3)*intmm(3))
791      Econstr=Econstr-magcon_lambda*(mag_1atom(1)*mag_1atom(1)+&
792 &     mag_1atom(2)*mag_1atom(2)+&
793 &     mag_1atom(3)*mag_1atom(3)) &
794 &     +magcon_lambda*(spinat(1,iatom)*spinat(1,iatom)+&
795 &     spinat(2,iatom)*spinat(2,iatom)+&
796 &     spinat(3,iatom)*spinat(3,iatom))
797      lVp=lVp+2*magcon_lambda*(intmm(1)*mag_1atom(1)+intmm(2)*mag_1atom(2)+intmm(3)*mag_1atom(3))
798    end if
799 
800    write(message, *) 'atom             constraining magnetic field'
801    call wrtout(std_out,message,'COLL')
802    write(message, '(I3,A2,E12.5,A2,E12.5,A2,E12.5)') &
803    iatom,'  ',magcon_lambda*intmm(1),'  ',magcon_lambda*intmm(2),'  ',magcon_lambda*intmm(3)
804    call wrtout(std_out,message,'COLL')
805 
806 !  End loop over atoms
807 !  -------------------------------------------
808  end do
809 
810 !Printing
811  write(message, '(A17,E10.3)' ) ' magcon_lambda    = ',magcon_lambda
812  call wrtout(std_out,message,'COLL')
813  write(message, '(A17,E12.5)' ) ' Lagrange penalty = ',Epen
814  call wrtout(std_out,message,'COLL')
815  write(message, '(A17,E12.5)' ) ' E_constraint     = ',Econstr
816  call wrtout(std_out,message,'COLL')
817  write(message, '(A17,E12.5)' ) ' lVp = ',lVp
818  call wrtout(std_out,message,'COLL')
819 
820  ABI_DEALLOCATE (intgden)
821 
822 end subroutine mag_constr_e

m_dens/radsmear [ Functions ]

[ Top ] [ m_dens ] [ Functions ]

NAME

 radsmear

FUNCTION

INPUTS

OUTPUT

PARENTS

CHILDREN

SOURCE

1365 function radsmear(r, rsph, rsm)
1366 
1367 
1368 !This section has been created automatically by the script Abilint (TD).
1369 !Do not modify the following lines by hand.
1370 #undef ABI_FUNC
1371 #define ABI_FUNC 'radsmear'
1372 !End of the abilint section
1373 
1374  implicit none
1375 
1376 !Arguments ------------------------------------
1377 !scalars
1378  real(dp) :: radsmear
1379  real(dp), intent(in) :: r, rsph, rsm
1380 
1381 !Local variables ------------------------------
1382 !scalars
1383  real(dp) :: xx
1384 
1385 !******************************************************************
1386 
1387  radsmear = zero
1388  if (r < rsph - rsm - tol12) then
1389    radsmear = one
1390  else if (r < rsph - tol12) then
1391    xx = (rsph - r) / rsm
1392    radsmear = xx**2*(3+xx*(1+xx*(-6+3*xx)))
1393  end if
1394 
1395 end function radsmear