TABLE OF CONTENTS
ABINIT/aim_follow [ Functions ]
NAME
aim_follow
FUNCTION
This routine follows the gradient line starting from the point vv. It stop when it arrives to the atom (nearer than rminl(iat)) or - if srch=true - also if it arrives under the already known part of Bader surface
COPYRIGHT
Copyright (C) 2002-2017 ABINIT group (PCasek,FF,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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt . WARNING This file does not follow the ABINIT coding rules (yet)
INPUTS
aim_dtset= the structured entity containing all input variables iatinit,iposinit= indexes of initial atom npmax= maximum number of division in each step
OUTPUT
iat,ipos= index of final atom nstep= returns the number of step needed
SIDE EFFECTS
srch= (true/false) check if the line is outside or
inside the atomic surface.
vv(3)= initial point in orthogonal coordinates
PARENTS
cpdrv,drvaim,rsurf
CHILDREN
critic,onestep,timein,vgh_rho
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 51 subroutine aim_follow(aim_dtset,vv,npmax,srch,iatinit,iposinit,iat,ipos,nstep) 52 53 use defs_basis 54 use defs_aimprom 55 use defs_parameters 56 use defs_abitypes 57 use m_errors 58 use m_profiling_abi 59 60 !This section has been created automatically by the script Abilint (TD). 61 !Do not modify the following lines by hand. 62 #undef ABI_FUNC 63 #define ABI_FUNC 'aim_follow' 64 use interfaces_18_timing 65 use interfaces_63_bader, except_this_one => aim_follow 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 integer,intent(in) :: iatinit,iposinit,npmax 73 integer,intent(out) :: iat,ipos,nstep 74 logical,intent(inout) :: srch 75 type(aim_dataset_type),intent(in) :: aim_dtset 76 !arrays 77 real(dp),intent(inout) :: vv(3) 78 79 !Local variables ------------------------------ 80 !scalars 81 integer :: i1,i2,i3,ii,iph,ires,ith,jj,kk,nit,np,nph,nsi,nth 82 real(dp) :: deltar,dg,dist,dph,dth,fac2,facf,h0old,hh,hold,rho,rr,rsmed 83 real(dp) :: t1,t2,t3,vcth,vph,vth,wall,xy,xyz 84 logical :: fin,ldebold,srchold,stemp,stemp2 85 character(len=50) :: formpc 86 character(len=500) :: msg 87 !arrays 88 real(dp) :: ev(3),grho(3),hrho(3,3),pom(3),vold(3),vt(3),vt1(3) 89 real(dp) :: zz(3,3) 90 91 !************************************************************************ 92 formpc='(":CP",2I5,3F12.8,3E12.4,I4,2E12.4)' 93 94 95 fin=.false. 96 97 srchold=srch 98 ldebold=ldeb 99 h0old=h0 100 101 nth=aim_dtset%nth 102 nph=aim_dtset%nph 103 104 if (slc==0) then 105 rminl(:)=aim_dtset%rmin 106 end if 107 108 if (deb) then 109 ldeb=.true. 110 end if 111 112 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 113 114 !Initial tests 115 116 if (iat/=0) then 117 if (rr<rminl(iat)) then 118 fin=.true. 119 write(std_out,*) 'rr < rmin iat=',iat,' ipos=',ipos 120 elseif (rho<aim_rhomin) then 121 fin=.true. 122 write(std_out,*) 'CHARGE LT rhomin ',rho,' < ',aim_rhomin 123 if (rho<zero) then 124 MSG_ERROR('RHO < 0 !!!') 125 end if 126 end if 127 end if 128 129 facf=aim_fac0 130 hh=aim_hmax 131 132 call timein(t1,wall) 133 nstep=0 134 nsi=0 135 136 !the principal cycle 137 138 madw : do while(.not.fin) 139 hold=hh 140 141 dg=vnorm(grho,0) 142 if (ldeb.or.deb) write(std_out,*) 'dg= ',dg 143 144 ! the time test 145 146 call timein(t3,wall) 147 t2=t3-t1 148 if (t2>300.0) then 149 write(std_out,*) 'TIME EXCEEDED 5 min IN FOLLOW' 150 write(std_out,*) 'h0 =',h0,' h =',hh,' h0old =',h0old,' dg =',dg 151 write(std_out,*) 'facf =',facf 152 msg = 'TIME EXCEEDED 5 min IN FOLLOW' 153 MSG_ERROR(msg) 154 end if 155 156 if (dg<aim_dgmin) then 157 write(std_out,*) 'gradient < dgmin ',dg,' < ',aim_dgmin 158 fin=.true. 159 iat=0 160 ipos=0 161 ! testing for the CP 162 if (npc>0) then 163 call critic(aim_dtset,vv,ev,zz,aim_dmaxcrit,ires,0) 164 if (ires==0) then 165 do jj=1,npc 166 pom(:)=pc(:,jj)-vv(:)+xatm(:,aim_dtset%batom) 167 dist=vnorm(pom,0) 168 if (dist<aim_tiny) cycle madw 169 end do 170 write(std_out,*) 'C.P. found !!' 171 npc=npc+1 172 do jj=1,3 173 pc(jj,npc)=vv(jj) 174 evpc(jj,npc)=ev(jj) 175 do kk=1,3 176 zpc(kk,jj,npc)=zz(kk,jj) 177 end do 178 end do 179 i1=ev(1)/abs(ev(1)) 180 i2=ev(2)/abs(ev(2)) 181 i3=ev(3)/abs(ev(3)) 182 icpc(npc)=i1+i2+i3 183 if (icpc(npc)==-3) then ! pseudoatom handling 184 npcm3=npcm3+1 185 write(std_out,*) 'Pseudo-atom found !!' 186 end if 187 188 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 189 write(22,formpc) 0,0,(pcrb(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),& 190 & ev(1)+ev(2)+ev(3),rho 191 write(std_out,formpc) 0,0,(pcrb(jj,npc),jj=1,3),(ev(jj),jj=1,3),icpc(npc),& 192 & ev(1)+ev(2)+ev(3),rho 193 else 194 write(std_out,*) 'C.P. not found !!' 195 end if 196 end if 197 198 cycle madw 199 end if 200 201 hh=h0/dg 202 if (ldeb.or.deb) write(std_out,*) 'h= ',hh,' h0= ',h0,' dg= ',dg 203 if (hh>aim_hmax) hh=aim_hmax 204 ! step modifications 205 206 hh=hh*facf 207 if (hh>(hold*aim_hmult)) then 208 hh=hold*aim_hmult 209 end if 210 211 do ii=1,3 212 vold(ii)=vv(ii) 213 end do 214 215 nit=0 216 hold=hh 217 218 ! one step following the gradient line 219 ! 220 call onestep(vv,rho,grho,hh,np,npmax,deltar) 221 do while (((np>npmax).or.(deltar>aim_stmax)).and.(deltar>aim_dmin)) 222 nit=nit+1 223 if (nit>5) then 224 if (deltar>aim_stmax) then 225 write(std_out,*) 'nit > 5 and deltar > stmax nit=',nit 226 else 227 write(std_out,*) 'nit > 5 and np > npmax nit=',nit 228 end if 229 end if 230 do ii=1,3 231 vv(ii)=vold(ii) 232 end do 233 hh=hh*0.3 234 call onestep(vv,rho,grho,hh,np,npmax,deltar) 235 end do 236 237 238 nstep=nstep+1 239 if (ldeb.or.deb) write(std_out,*) 'h= ',hh 240 241 fac2=hh/hold 242 if (fac2>=1._dp) then 243 facf=facf*1.2 244 else 245 if (fac2>=aim_facmin) then 246 facf=fac2 247 else 248 facf=aim_facmin 249 end if 250 end if 251 252 if (deb.or.ldeb) then 253 write(std_out,*) ':POS ',vv 254 write(std_out,*) ':RBPOS ',vt1 255 write(std_out,*) ':GRAD ',grho 256 end if 257 258 call vgh_rho(vv,rho,grho,hrho,rr,iat,ipos,slc) 259 dg=vnorm(grho,0) 260 pom(:)=vv(:)-xatm(:,iatinit)-atp(:,iposinit) 261 262 if (iat /= 0) then 263 fin=.true. 264 write(std_out,*) 'r < rmin iat=',iat,' ipos=',ipos 265 cycle madw 266 end if 267 268 if (rho<aim_rhomin) then 269 fin=.true. 270 write(std_out,*) 'charge < rhomin ',rho,' < ',aim_rhomin 271 if (rho<zero) then 272 MSG_ERROR('RHO < 0 !!!') 273 end if 274 iat=0 275 ipos=0 276 cycle madw 277 end if 278 279 if (npcm3>0) then 280 do jj=1,npc 281 if (icpc(jj)==(-3)) then 282 pom(:)=pc(:,jj)-vv(:)+xatm(:,aim_dtset%batom) 283 dist=vnorm(pom,0) 284 if (dist<(aim_dtset%rmin**2*0.1)) then 285 iat=0 286 ipos=0 287 fin=.true. 288 write(std_out,*) 'We are inside a pseudo-atom' 289 cycle madw 290 end if 291 end if 292 end do 293 end if 294 295 nsi=nsi+1 296 297 ! surface checking 298 299 if (srch.and.(nsi>=nsimax)) then 300 nsi=0 301 ith=0 302 iph=0 303 do ii=1,3 304 vt(ii)=vv(ii)-xatm(ii,iatinit) 305 end do 306 xy=vt(1)*vt(1)+vt(2)*vt(2) 307 xyz=xy+vt(3)*vt(3) 308 xyz=sqrt(xyz) 309 if (xy<aim_snull) then 310 vcth=1._dp 311 if (vt(3)<0._dp) vcth=-vcth 312 vph=0._dp 313 else 314 vcth=vt(3)/xyz 315 vph=atan2(vt(2),vt(1)) 316 end if 317 vth=acos(vcth) 318 if (vth<th(1)) then 319 ith=0 320 else 321 if (vth>th(nth)) then 322 ith=nth 323 else 324 do ii=2,nth 325 if (vth<th(ii)) then 326 ith=ii-1 327 exit 328 end if 329 end do 330 end if 331 end if 332 333 if (vph<ph(1)) then 334 iph=0 335 else 336 if (vph>ph(nph)) then 337 iph=nph 338 else 339 do ii=2,nph 340 if (vph<ph(ii)) then 341 iph=ii-1 342 exit 343 end if 344 end do 345 end if 346 end if 347 348 stemp=(iph>0).and.(iph<nph) 349 stemp=stemp.and.(ith>0).and.(ith<nth) 350 351 if (stemp) then 352 stemp2=rs(ith,iph)>0._dp 353 stemp2=stemp2.and.(rs(ith+1,iph)>0._dp) 354 stemp2=stemp2.and.(rs(ith+1,iph+1)>0._dp) 355 stemp2=stemp2.and.(rs(ith,iph+1)>0._dp) 356 if (stemp2) then 357 dth=th(ith+1)-th(ith) 358 dph=ph(iph+1)-ph(iph) 359 rsmed=rs(ith,iph)*(th(ith+1)-vth)/dth*(ph(iph+1)-vph)/dph 360 rsmed=rsmed+rs(ith+1,iph)*(vth-th(ith))/dth*(ph(iph+1)-vph)/dph 361 rsmed=rsmed+rs(ith+1,iph+1)*(vth-th(ith))/dth*(vph-ph(iph))/dph 362 rsmed=rsmed+rs(ith,iph+1)*(th(ith+1)-vth)/dth*(vph-ph(iph))/dph 363 if (rsmed>xyz) then 364 write(std_out,*) 'We are inside the surface' 365 iat=iatinit 366 ipos=iposinit 367 else 368 write(std_out,*) 'We are outside the surface' 369 iat=0 370 ipos=0 371 end if 372 fin=.true. 373 cycle madw 374 end if 375 end if 376 end if 377 378 end do madw 379 380 381 srch=srchold 382 ldeb=ldebold 383 h0=h0old 384 385 386 end subroutine aim_follow