./67_common/spatialchempot.F90

ABINIT/spatialchempot

ABINIT/spatialchempot [ Modules ]

NAME

  spatialchempot

FUNCTION

  Treat spatially varying chemical potential.
  Compute energy and derivative with respect to dimensionless reduced atom coordinates of the
  spatially varying chemical potential. No contribution to stresses.

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

INPUTS

 chempot(3,nzchempot,ntypat)=input array with information about the chemical potential (see input variable description)
 natom=number of atoms in unit cell
 ntypat=number of type of atoms
 nzchempot=number of limiting planes for chemical potential
 typat(natom)=integer label of each type of atom (1,2,...)
 xred(3,natom)=relative coords of atoms in unit cell (dimensionless)

OUTPUT

 e_chempot=chemical potential energy in hartrees
 grchempottn(3,natom)=grads of e_chempot wrt xred(3,natom), hartrees.

PARENTS

      setvtr

CHILDREN

SOURCE

 36 #if defined HAVE_CONFIG_H
 37 #include "config.h"
 38 #endif
 39 
 40 #include "abi_common.h"
 41 
 42  subroutine spatialchempot(e_chempot,chempot,grchempottn,natom,ntypat,nzchempot,typat,xred)
 43 
 44  use defs_basis
 45 
 46 !This section has been created automatically by the script Abilint (TD).
 47 !Do not modify the following lines by hand.
 48 #undef ABI_FUNC
 49 #define ABI_FUNC 'spatialchempot'
 50 !End of the abilint section
 51 
 52  implicit none
 53 
 54 !Arguments ------------------------------------
 55 !scalars
 56  integer,intent(in) :: natom,ntypat,nzchempot
 57  real(dp),intent(out) :: e_chempot
 58 !arrays
 59  integer,intent(in) :: typat(natom)
 60  real(dp),intent(in) :: chempot(3,nzchempot,ntypat),xred(3,natom)
 61  real(dp),intent(out) :: grchempottn(3,natom)
 62 
 63 !Local variables-------------------------------
 64 !scalars
 65  integer :: iatom,itypat,iz
 66  real(dp) :: a_2,a_3,cp0,cp1,dcp0,dcp1,ddz,deltaz,deltaziz
 67  real(dp) :: dqz,dz1,qz,zred,z0
 68 !character(len=500) :: message
 69 
 70 ! *************************************************************************
 71 
 72 !DEBUG
 73 !write(std_out,'(a)')' spatialchempot : enter '
 74 !write(std_out,'(a,2i6)')' nzchempot,ntypat',nzchempot,ntypat
 75 !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,1)=',chempot(1:3,1:2,1)
 76 !write(std_out,'(a,6es13.3)') ' chempot(1:3,1:2,2)=',chempot(1:3,1:2,2)
 77 !ENDDEBUG
 78 
 79  e_chempot=zero
 80  grchempottn(:,:)=zero
 81 
 82 !Loop on the different atoms
 83  do iatom=1,natom
 84 
 85    itypat=typat(iatom)
 86    zred=xred(3,iatom)
 87 
 88 !  Determine the delimiting plane just lower to zred
 89 !  First compute the relative zred with respect to the first delimiting plane
 90 !  Take into account a tolerance : 
 91    deltaz=zred-chempot(1,1,itypat)
 92    deltaz=modulo(deltaz+tol12,1.0d0)-tol12
 93 !  deltaz is positive (or higher than -tol12), and lower than one-tol12.
 94    do iz=2,nzchempot+1
 95      if(iz/=nzchempot+1)then
 96        deltaziz=chempot(1,iz,itypat)-chempot(1,1,itypat)
 97      else
 98        deltaziz=one
 99      end if
100      if(deltaziz>deltaz)exit
101    end do
102 
103 !  Defines coordinates and values inside the delimiting interval, 
104 !  with respect to the lower delimiting plane
105    z0=chempot(1,iz-1,itypat)-chempot(1,1,itypat) ; cp0=chempot(2,iz-1,itypat) ; dcp0=chempot(3,iz-1,itypat)
106    if(iz/=nzchempot+1)then
107      dz1=chempot(1,iz,itypat)-chempot(1,iz-1,itypat) ; cp1=chempot(2,iz,itypat) ; dcp1=chempot(3,iz,itypat)
108    else
109      dz1=(chempot(1,1,itypat)+one)-chempot(1,nzchempot,itypat) ; cp1=chempot(2,1,itypat) ; dcp1=chempot(3,1,itypat)
110    end if
111    ddz=deltaz-z0
112 
113 !DEBUG
114 !  write(std_out,'(a,2i5)')' Delimiting planes, iz-1 and iz=',iz-1,iz
115 !  write(std_out,'(a,2es13.3)')' z0,  dz1= :',z0,dz1
116 !  write(std_out,'(a,2es13.3)')' cp0, cp1= :',cp0,cp1
117 !  write(std_out,'(a,2es13.3)')' dcp0, dcp1= :',dcp0,dcp1
118 !  write(std_out,'(a,2es13.3)')' deltaz,ddz=',deltaz,ddz
119 !ENDDEBUG
120 
121 !  Determine the coefficient of the third-order polynomial taking z0 as origin
122 !  P(dz=z-z0)= a_3*dz**3 + a_2*dz**2 + a_1*dz + a_0 ; obviously a_0=cp0 and a_1=dcp0
123 !  Define qz=a_3*dz + a_2 and dqz=3*a_3*dz + 2*a_2
124    qz=((cp1-cp0)-dcp0*dz1)/dz1**2
125    dqz=(dcp1-dcp0)/dz1
126    a_3=(dqz-two*qz)/dz1
127    a_2=three*qz-dqz
128 
129 !  Compute value and gradient of the chemical potential, at ddz wrt to lower delimiting plane
130    e_chempot=e_chempot+(a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0)
131    grchempottn(3,iatom)=three*a_3*ddz**2 + two*a_2*ddz + dcp0
132 
133 !DEBUG
134 !  write(std_out,'(a,4es16.6)')' qz,dqz=',qz,dqz
135 !  write(std_out,'(a,4es16.6)')' cp0,dcp0,a_2,a_3=',cp0,dcp0,a_2,a_3
136 !  write(std_out,'(a,2es13.3)')' dcp0*ddz + cp0=',dcp0*ddz + cp0
137 !  write(std_out,'(a,2es13.3)')' a_2*ddz**2=',a_2*ddz**2
138 !  write(std_out,'(a,2es13.3)')' a_3*ddz**3=',a_3*ddz**3
139 !  write(std_out,'(a,2es13.3)')' contrib=',a_3*ddz**3 + a_2*ddz**2 + dcp0*ddz + cp0
140 !  write(std_out,'(a,2es13.3)')' e_chempot=',e_chempot
141 !  write(std_out,'(a,3es20.10)')' grchempottn=',grchempottn(:,iatom)
142 !ENDDEBUG
143 
144  end do
145 
146 !DEBUG
147 !write(std_out,'(a)')' spatialchempot : exit '
148 !write(std_out,'(a,es16.6)') ' e_chempot=',e_chempot
149 !ENDDEBUG
150 
151 end subroutine spatialchempot

TABLE OF CONTENTS

ABINIT/spatialchempot [ Modules ]