TABLE OF CONTENTS
ABINIT/chkpawovlp [ Functions ]
NAME
chkpawovlp
FUNCTION
Verify that the paw spheres are not overlapping
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (FJ, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
natom=number of atoms in cell.
ntypat=number of types of atoms in unit cell.
pawovlp=percentage of voluminal overlap ratio allowed to continue execution
(if negative value, execution always continues)
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
typat(natom)=type (integer) for each atom
xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
(only checking)
NOTES
PARENTS
bethe_salpeter,respfn,scfcv,screening,sigma,wfk_analyze
CHILDREN
wrtout
SOURCE
39 #if defined HAVE_CONFIG_H 40 #include "config.h" 41 #endif 42 43 #include "abi_common.h" 44 45 subroutine chkpawovlp(natom,ntypat,pawovlp,pawtab,rmet,typat,xred) 46 47 use defs_basis 48 use m_errors 49 use m_profiling_abi 50 51 use m_pawtab, only : pawtab_type 52 53 !This section has been created automatically by the script Abilint (TD). 54 !Do not modify the following lines by hand. 55 #undef ABI_FUNC 56 #define ABI_FUNC 'chkpawovlp' 57 use interfaces_14_hidewrite 58 !End of the abilint section 59 60 implicit none 61 62 !Arguments --------------------------------------------- 63 !scalars 64 integer,intent(in) :: natom,ntypat 65 real(dp) :: pawovlp 66 !arrays 67 integer,intent(in) :: typat(natom) 68 real(dp),intent(in) :: rmet(3,3),xred(3,natom) 69 type(pawtab_type),intent(in) :: pawtab(ntypat) 70 71 !Local variables --------------------------------------- 72 !scalars 73 integer :: ia,ib,ii,t1,t2,t3 74 logical :: stop_on_error 75 real(dp) :: dd,dif1,dif2,dif3,ha,hb,norm2 76 real(dp) :: ratio_percent,va,vb,vv 77 character(len=750) :: message 78 !arrays 79 integer :: iamax(2),ibmax(2),iovl(2) 80 real(dp) :: norm2_min(2),r2cut(2),ratio_percent_max(2),rcuta(2),rcutb(2) 81 82 83 ! ************************************************************************* 84 85 DBG_ENTER("COLL") 86 87 iamax(:)=-1;ibmax(:)=-1 88 norm2_min(:)=-1.d0;ratio_percent_max(:)=-1.d0 89 iovl(:)=0 90 91 !Loop on "overlapping" atoms with the maximum overlap 92 do ia=1,natom 93 94 rcuta(1)=pawtab(typat(ia))%rpaw 95 rcuta(2)=pawtab(typat(ia))%rshp 96 97 do ib=ia,natom 98 99 rcutb(1)=pawtab(typat(ib))%rpaw 100 rcutb(2)=pawtab(typat(ib))%rshp 101 r2cut(1)=(rcuta(1)+rcutb(1))**2 102 r2cut(2)=(rcuta(2)+rcutb(2))**2 103 104 ! Visit the box and its first images: 105 do t3=-1,1 106 do t2=-1,1 107 do t1=-1,1 108 109 dif1=xred(1,ia)-(xred(1,ib)+dble(t1)) 110 dif2=xred(2,ia)-(xred(2,ib)+dble(t2)) 111 dif3=xred(3,ia)-(xred(3,ib)+dble(t3)) 112 norm2=sqnrm_pawovlp(dif1,dif2,dif3) 113 114 do ii=1,2 115 116 if(norm2>tol10.and.norm2<r2cut(ii)) then 117 118 iovl(ii)=iovl(ii)+1 119 120 ! Compute the overlap ratio: 121 dd=sqrt(norm2) 122 va=4._dp/3._dp*pi*rcuta(ii)**3 123 vb=4._dp/3._dp*pi*rcutb(ii)**3 124 ha=(rcutb(ii)**2-(dd-rcuta(ii))**2)/(two*dd) 125 hb=(rcuta(ii)**2-(dd-rcutb(ii))**2)/(two*dd) 126 vv=pi/3.d0*(ha**2*(three*rcuta(ii)-ha)+hb**2*(three*rcutb(ii)-hb)) 127 ratio_percent=100._dp*min(vv/min(va,vb),one) 128 if (ratio_percent>ratio_percent_max(ii)) then 129 ratio_percent_max(ii)=ratio_percent 130 norm2_min(ii)=norm2 131 iamax(ii)=ia;ibmax(ii)=ib 132 end if 133 134 end if 135 end do 136 end do 137 end do 138 end do 139 end do 140 end do 141 142 stop_on_error=(abs(pawovlp)<=tol6.or.(pawovlp>tol6.and.ratio_percent_max(1)>pawovlp)) 143 144 !Print adapted message with overlap value 145 if (iovl(1)+iovl(2)>0) then 146 147 !ii=1: PAW augmentation regions overlap 148 !ii=2: compensation charges overlap 149 if (iovl(2)==0) ii=1 150 if (iovl(2)> 0) ii=2 151 152 if (iovl(ii)>0) then 153 154 if (ii==1) write(message,' (a)' ) 'PAW SPHERES ARE OVERLAPPING!' 155 if (ii==2) write(message, '(2a)' )'PAW COMPENSATION DENSITIES ARE OVERLAPPING !!!!' 156 157 if (iovl(ii)==1) then 158 write(message, '(3a)' ) trim(message),ch10,& 159 & ' There is one pair of overlapping atoms.' 160 else 161 write(message, '(3a,i5,a)' ) trim(message),ch10,& 162 & ' There are ', iovl(1),' pairs of overlapping atoms.' 163 end if 164 write(message, '(3a,i3,a,i3,a)' ) trim(message),ch10,& 165 ' The maximum overlap percentage is obtained for the atoms ',iamax(ii),' and ',ibmax(ii),'.' 166 write(message, '(2a,2(a,i3),a,f9.5,a,2(a,i3,a,f9.5,a),a,f5.2,a)' ) trim(message),ch10,& 167 & ' | Distance between atoms ',iamax(ii),' and ',ibmax(ii),' is : ',sqrt(norm2_min(ii)),ch10,& 168 & ' | PAW radius of the sphere around atom ',iamax(ii),' is: ',pawtab(typat(iamax(ii)))%rpaw,ch10,& 169 & ' | PAW radius of the sphere around atom ',ibmax(ii),' is: ',pawtab(typat(ibmax(ii)))%rpaw,ch10,& 170 & ' | This leads to a (voluminal) overlap ratio of ',ratio_percent_max(ii),' %' 171 if (ii==2) then 172 write(message, '(3a)' ) trim(message),ch10,& 173 & 'THIS IS DANGEROUS !, as PAW formalism assumes non-overlapping compensation densities.' 174 end if 175 176 if (stop_on_error) then 177 MSG_ERROR_NOSTOP(message,ia) !ia is dummy 178 else 179 MSG_WARNING(message) 180 end if 181 182 end if 183 184 ! Print advice 185 if (stop_on_error) then 186 write(message, '(3a)' )& 187 & ' Action: 1- decrease cutoff radius of PAW dataset',ch10,& 188 & ' OR 2- ajust "pawovlp" input variable to allow overlap (risky)' 189 MSG_ERROR(message) 190 end if 191 192 ! Print last message if execution continues: 193 if (pawovlp<=tol6) then 194 write(message, '(6a)' ) & 195 & ' Results might be approximate,',ch10,& 196 & ' and even inaccurate (if overlap is too big) !',ch10,& 197 & ' Assume experienced user. Execution will continue.',ch10 198 call wrtout(std_out,message,'COLL') 199 else if (ratio_percent_max(1)<=pawovlp) then 200 write(message, '(8a)' ) & 201 & ' Overlap ratio seems to be acceptable (less than value',ch10,& 202 & ' of "pawovlp" input parameter): execution will continue.',ch10,& 203 & ' But be aware that results might be approximate,',ch10,& 204 & ' and even inaccurate (depending on your physical system) !',ch10 205 call wrtout(std_out,message,'COLL') 206 end if 207 208 end if !iovl>0 209 210 DBG_EXIT("COLL") 211 212 contains 213 214 function sqnrm_pawovlp(u1,u2,u3) 215 !squared norm of a vector 216 217 !This section has been created automatically by the script Abilint (TD). 218 !Do not modify the following lines by hand. 219 #undef ABI_FUNC 220 #define ABI_FUNC 'sqnrm_pawovlp' 221 !End of the abilint section 222 223 real(dp) :: sqnrm_pawovlp 224 real(dp),intent(in) :: u1,u2,u3 225 226 sqnrm_pawovlp=rmet(1,1)*u1*u1+rmet(2,1)*u2*u1+rmet(3,1)*u3*u1& 227 & +rmet(1,2)*u1*u2+rmet(2,2)*u2*u2+rmet(3,2)*u3*u2& 228 & +rmet(1,3)*u1*u3+rmet(2,3)*u2*u3+rmet(3,3)*u3*u3 229 230 end function sqnrm_pawovlp 231 end subroutine chkpawovlp