TABLE OF CONTENTS

ABINIT/rsurf [ Functions ]

[ Top ] [ Functions ]

NAME

 rsurf

FUNCTION

 Basic routine for determination of the radius of Bader surface
 for spherical rayon theta,phi
 the bassin is tested by following the gradient line
 If srch==true (in general for calls from surf) the routine aim_follow
 is called to stop when it arrives under already known part of surface
 Simple bissection method is used to obtain the radius

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
 rr0= starting radius
 theta,phi = the spherical direction
 iatinit= the atom index
 srch= see above
 npmax= maximum number of divisions in one step for follow

OUTPUT

 rr= radius
 grho(3)= gradient on the surface

PARENTS

      drvaim,surf

CHILDREN

      aim_follow,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 rsurf(aim_dtset,rr,grho,theta,phi,rr0,iatinit,npmax,srch)
 52 
 53  use m_profiling_abi
 54 
 55  use defs_basis
 56  use defs_parameters
 57  use defs_aimprom
 58  use defs_abitypes
 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 'rsurf'
 64  use interfaces_18_timing
 65  use interfaces_63_bader, except_this_one => rsurf
 66 !End of the abilint section
 67 
 68  implicit none
 69 
 70 !Arguments ------------------------------------
 71 !scalars
 72  integer,intent(in) :: iatinit,npmax
 73  real(dp),intent(in) :: phi,rr0,theta
 74  real(dp),intent(out) :: rr
 75  logical,intent(in) :: srch
 76 !arrays
 77  real(dp),intent(out) :: grho(3)
 78 !no_abirules
 79  type(aim_dataset_type),intent(in) :: aim_dtset
 80 
 81 !Local variables ------------------------------
 82 !scalars
 83  integer :: iat,ii,ipos,iposinit,jj,nstep
 84  real(dp),parameter :: mfkt=1.d1
 85  real(dp) :: aa,dmax,dr,drr,rho,rr1,rr2,t1,t2,wall
 86  logical :: cross,deb_tmp,in,in1,in2,low,srch_tmp
 87 !arrays
 88  real(dp) :: hrho(3,3),unvec(3),vv(3)
 89 
 90 ! *********************************************************************
 91 
 92  srch_tmp=srch
 93  deb_tmp=deb
 94 
 95 !unity vecteur in the direction (theta,phi)
 96 
 97  unvec(1)=sin(theta)*cos(phi)
 98  unvec(2)=sin(theta)*sin(phi)
 99  unvec(3)=cos(theta)
100 
101 
102  rr=rr0
103  rr1=rr
104  rr2=rr
105  drr=1._dp
106  if (abs(rr0-r0)<1.0d-12) then
107    dr=aim_dtset%dr0*mfkt
108  else
109    dr=aim_dtset%dr0
110  end if
111 
112  vv(1)=xatm(1,aim_dtset%batom)
113  vv(2)=xatm(2,aim_dtset%batom)
114  vv(3)=xatm(3,aim_dtset%batom)
115 
116 
117  iposinit=batcell
118  write(std_out,'("ATOM iat=",i4," ipos=",i4)') aim_dtset%batom,batcell
119  jj=0
120 
121  cross=.false.
122 
123  in=.true.
124  low=.false.
125 
126  dmax=h0
127 
128  in1=.true.
129  in2=in1
130 
131  do while((drr>aim_drmin).or.(jj<2))
132    call timein(t1,wall)
133    jj=jj+1
134    do ii=1,3
135      vv(ii)=xatm(ii,aim_dtset%batom)+rr*unvec(ii)
136    end do
137 
138 !  VACUUM CONDITION
139 
140    call vgh_rho(vv,rho,grho,hrho,aa,iat,ipos,0)
141    if (rho < aim_rhomin) exit
142 
143    ldeb=.false.
144 
145    call aim_follow(aim_dtset,vv,npmax,srch_tmp,iatinit,iposinit,iat,ipos,nstep)
146 
147    call timein(t2,wall)
148    t2=t2-t1
149 
150    write(std_out,'(a,i4,a,f12.8,a,i4,a,i4,a,f10.5,a,i4)') &
151 &   ' :STEP ',jj,' r=',rr,' iat=',iat,' ipos=',ipos,' time(sec)=',t2,' nstep=',nstep
152 
153    if ((iat.eq.iatinit).and.(ipos.eq.iposinit)) then
154      in=.true.
155    else
156      in=.false.
157    end if
158 
159 !  
160 !  NEW RADIUS
161 !  
162 
163    if ((jj.eq.1).or.((in1.eqv.in).and.(.not.cross))) then
164      if (in) then
165        rr2=rr1
166        rr1=rr
167        rr=rr+dr
168      else
169        rr2=rr1
170        rr1=rr
171        rr=rr-dr
172      end if
173      if ((jj>2).and.(dr<(0.6))) then
174 !      modification of the step
175        dr=dr*aim_fac
176        if (deb_tmp) write(std_out,*) ':DR ',dr
177      end if
178    else
179      if (.not.cross) then
180        cross=.true.
181        rr2=rr1
182      else
183        if (in2) then
184          if (in) then
185            rr2=rr1
186          else
187            in1=in2
188          end if
189        else
190          if (in) then
191            in1=in2
192          else
193            rr2=rr1
194          end if
195        end if
196      end if
197      rr1=rr
198      rr=(rr2+rr1)/2.0
199    end if
200 
201    in2=in1
202    in1=in
203    drr=abs(rr2-rr1)/rr
204    if (deb_tmp) write(std_out,*) ':DRR ',jj,rr2,rr1,drr
205  end do
206 
207 end subroutine rsurf