TABLE OF CONTENTS


ABINIT/critic [ Functions ]

[ Top ] [ Functions ]

NAME

 critic

FUNCTION

     Search for a critical point starting from point vv

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
 dmax= maximal step
 sort= 0(default) general CP searching (Newton-Raphson)
                  -1,1,3 searching of specific type CP (Popelier)

OUTPUT

 ev= eigenvalues (ordered) of the Hessian in the final point
 zz=  eigenvectors of the Hessian in the final point
 ires= if ires==0 => CP found
       if ires==1 => CP not found within the maximum steps

SIDE EFFECTS

 vv(3)= starting point and final point

PARENTS

      aim_follow,cpdrv,critics

CHILDREN

SOURCE

 41 #if defined HAVE_CONFIG_H
 42 #include "config.h"
 43 #endif
 44 
 45 #include "abi_common.h"
 46 
 47 subroutine critic(aim_dtset,vv,ev,zz,dmax,ires,sort)
 48 
 49  use defs_basis
 50  use defs_aimprom
 51  use defs_parameters
 52  use defs_abitypes
 53  use m_errors
 54  use m_profiling_abi
 55 
 56 !This section has been created automatically by the script Abilint (TD).
 57 !Do not modify the following lines by hand.
 58 #undef ABI_FUNC
 59 #define ABI_FUNC 'critic'
 60  use interfaces_28_numeric_noabirule
 61  use interfaces_63_bader, except_this_one => critic
 62 !End of the abilint section
 63 
 64  implicit none
 65 
 66 !Arguments ------------------------------------
 67 !scalars
 68  integer,intent(in) :: sort
 69  integer,intent(out) :: ires
 70  real(dp),intent(in) :: dmax
 71 !arrays
 72  real(dp),intent(inout) :: vv(3)
 73  real(dp),intent(out) :: ev(3),zz(3,3)
 74 !no_abirules
 75  type(aim_dataset_type), intent(in) :: aim_dtset
 76 
 77 !Local variables ------------------------------
 78 !scalars
 79  integer :: iat,id,ii,info,ipos,istep,jii,jj,nrot
 80  real(dp),parameter :: evol=1.d-3
 81  real(dp) :: chg,dg,dltcmax,dv,dvold,rr,ss
 82  logical :: oscl,outof
 83 !arrays
 84  integer :: ipiv(3)
 85  real(dp) :: dc(3),ff(3),grho(3),hrho(3,3),lp(3),vold(3),vt(3),yy(3,3)
 86  real(dp),allocatable :: lamb(:),pom(:,:),pom2(:,:)
 87 
 88 !************************************************************************
 89 
 90 !DEBUG
 91 !write(std_out,*)' critic : enter '
 92 !ENDDEBUG
 93  oscl=.false.
 94  if (sort==3) then
 95    ABI_ALLOCATE(pom,(4,4))
 96    ABI_ALLOCATE(pom2,(4,4))
 97    ABI_ALLOCATE(lamb,(4))
 98  elseif (sort/=0) then
 99    ABI_ALLOCATE(pom,(3,3))
100    ABI_ALLOCATE(pom2,(3,3))
101    ABI_ALLOCATE(lamb,(3))
102  end if
103 
104 
105  deb=.false.
106  istep=0
107  ires=0
108 
109 !DEBUG
110 !write(std_out,'(":POSIN ",3F16.8)') vv
111 !do jj=1,3
112 !vt(jj)=rprimd(1,jj)*vv(1)+rprimd(2,jj)*vv(2)+rprimd(3,jj)*vv(3)
113 !end do
114 !write(std_out,'(":RBPOSIN ",3F16.8)') vt
115 !ENDDEBUG
116 
117  call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0)
118 
119 !write(std_out,'(":GRAD ",3F16.8)') grho
120 !write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jj),jj=1,3),ii=1,3)
121 
122  dg=1.0_dp
123  dv=1.0_dp
124  dg = vnorm(grho,0)
125 
126  if (chg < aim_rhomin) then
127    ires=1
128 !  DEBUG
129 !  write(std_out,*)' critic : exit, ires=1'
130 !  ENDDEBUG
131    return
132  end if
133 
134 !main cycle => limits (adhoc):
135 !aim_dtset%lstep - minimal step
136 !aim_dtset%lgrad - minimal norm of gradient
137 !aim_maxstep - max number of steps
138 
139  do while ((dv>aim_dtset%lstep).and.(dg>aim_dtset%lgrad).and.(istep<aim_maxstep))
140    istep=istep+1
141    vold(:)=vv(:)
142    dvold=dv
143    ev(:)=0._dp
144    yy(:,:)=0._dp
145    call jacobi(hrho,3,3,ev,yy,nrot)   ! eigenval of Hessian
146    call ordr(ev,yy,3,-1)  ! ordering
147 
148 !  modification of the Newton-Raphson step to searching
149 !  specific type of CP (Popelier algorithm)
150 
151    ff(:)=0._dp
152    lp(:)=0._dp
153    dc(:)=0._dp
154    outof=.false.
155    do ii=1,3
156      do jj=1,3
157        ff(ii)=ff(ii)+yy(jj,ii)*grho(jj)
158      end do
159    end do
160    id=sign(1._dp,ev(1))+sign(1._dp,ev(2))+sign(1._dp,ev(3))
161    if (id /= sort) then
162      outof=.true.
163      select case (sort)
164      case (-1)
165        lp(3)=0.5_dp*(ev(3)-sqrt(ev(3)*ev(3)+4._dp*ff(3)*ff(3)))
166        pom(:,:)=0._dp
167        pom2(:,:)=0._dp
168        lamb(:)=0._dp
169        do ii=1,2
170          pom(ii,ii)=ev(ii)
171          pom(ii,3)=ff(ii)
172          pom(3,ii)=ff(ii)
173        end do
174        call jacobi(pom,3,3,lamb,pom2,nrot)
175        call ordr(lamb,pom2,3,1)
176        do ii=1,3
177          lp(1)=lamb(ii)
178          if (abs(pom2(3,ii))>1.0d-24) exit
179        end do
180        lp(2)=lp(1)
181 
182 !        write(std_out,*) (ev(ii),ii=1,3)
183 !        write(std_out,*) (lamb(ii),ii=1,3)
184 !        write(std_out,*) ':ID  ',id,lp(1),lp(3)
185 
186      case (1)
187        lp(1)=0.5_dp*(ev(1)+sqrt(ev(1)*ev(1)+4._dp*ff(1)*ff(1)))
188        pom(:,:)=0._dp
189        pom2(:,:)=0._dp
190        lamb(:)=0._dp
191        do ii=2,3
192          pom(ii-1,ii-1)=ev(ii)
193          pom(ii-1,3)=ff(ii)
194          pom(3,ii-1)=ff(ii)
195        end do
196        call jacobi(pom,3,3,lamb,pom2,nrot)
197        call ordr(lamb,pom2,3,1)
198        do ii=3,1,-1
199          lp(2)=lamb(ii)
200          if (abs(pom2(3,ii))>1.0d-24) exit
201        end do
202        lp(3)=lp(2)
203 
204      case (3)
205        pom(:,:)=0._dp
206        pom2(:,:)=0._dp
207        lamb(:)=0._dp
208        do ii=1,3
209          pom(ii,ii)=ev(ii)
210          pom(ii,4)=ff(ii)
211          pom(4,ii)=ff(ii)
212        end do
213        call jacobi(pom,4,4,lamb,pom2,nrot)
214        call ordr(lamb,pom2,4,1)
215        do ii=4,1,-1
216          lp(1)=lamb(ii)
217          if (abs(pom2(4,ii))>1.0d-24) exit
218        end do
219        lp(2)=lp(1); lp(3)=lp(1)
220      case default
221        lp(:)=0._dp
222      end select
223    end if
224 
225    do ii=1,3
226      if (abs(ev(ii)-lp(ii))<1.0d-24) then
227        outof=.false.
228        exit
229      end if
230    end do
231    do ii=1,3                      ! SEARCHING STEP
232      do jj=1,3
233        if (outof) then
234          dc(ii)=dc(ii)+ff(jj)*yy(ii,jj)/(ev(jj)-lp(jj))
235        elseif (abs(ev(jj))>1.0d-24) then
236          dc(ii)=dc(ii)+ff(jj)*yy(ii,jj)/ev(jj)
237        else
238          MSG_ERROR("zero eigval of Hessian")
239        end if
240      end do
241    end do
242 
243    dltcmax = vnorm(dc,0)
244    if (dltcmax>dmax) then                 ! STEP RESTRICTION
245      do ii=1,3
246        dc(ii)=dc(ii)*dmax/dltcmax
247      end do
248    end if                                  ! primitive handling of oscillations
249    ss=vnorm(dc,0)                          ! usually not needed
250    ss=abs(ss-dv)/ss
251    if ((ss < evol).and.(oscl)) then
252      dc(:)=dc(:)/2._dp
253    end if
254 
255 
256    do ii=1,3
257      vv(ii) = vv(ii) - dc(ii)
258    end do
259 
260 !  DEBUG
261 !  write(std_out,'(":POSIN ",3F16.8)') vv
262 !  ENDDEBUG
263 
264    call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0)
265    dg = vnorm(grho,0)
266 
267    if (deb) then                 !  DEBUGG OUTPUT
268      write(std_out,'("AFTER STEP ===================================")')
269      write(std_out,'(":HESSIAN^(-1) ",/,3F16.8,/,3F16.8,/,3F16.8)') ((yy(ii,jii),jii=1,3),ii=1,3)
270      write(std_out,'(":DC ",3F16.8)') dc
271      write(std_out,*) 'STEP ',istep
272      write(std_out,'(":POS ",3F16.8)') vv
273      write(std_out,'(":GRAD ",3F16.8)') grho
274      write(std_out,*) ':DGRAD,CHG ',dg,chg
275      write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)') ((hrho(ii,jii),jii=1,3),ii=1,3)
276    end if
277    vt(:)=vv(:)-vold(:)
278    dv=vnorm(vt,0)
279    ss=abs(dvold-dv)/dv
280    if (ss < evol) oscl=.true.
281  end do
282 
283 !end of main cycle
284 
285 !the final output
286 
287  write(std_out,*) 'iste:',istep, dv, dg
288  if (istep>=aim_maxstep)then
289    write(std_out,*) ' istep=MAXSTEP ! Examine lstep2 and lgrad2 .'
290    if ( (dv>aim_dtset%lstep2) .and. (dg>aim_dtset%lgrad2 )) then
291      write(std_out,'(":POSOUT ",3F16.8)') vv
292      ires=1
293    end if
294  end if
295 
296  vt(:)=vv(:)
297 
298 !write(std_out,'(":POSOUT ",3F16.8)') vv
299 !write(std_out,'(":RBPOSOUT ",3F16.8)') vt
300 
301  call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0)
302 
303 !write(std_out,'(":GRAD ",3F16.8)') grho
304 !write(std_out,'(":HESSIAN ",/,3F16.8,/,3F16.8,/,3F16.8)')&
305 !& ((hrho(ii,jii),jii=1,3),ii=1,3)
306 
307 
308 !FINAL INVERSION OF HESSIAN
309 
310  call ludcmp(hrho,3,3,ipiv,id,info)
311  if (info /= 0) then
312    write(std_out,*) 'Error inverting hrho:'
313    do ii=1,3
314      write(std_out,*) (hrho(ii,jii),jii=1,3)
315    end do
316    ires=1
317 !  DEBUG
318 !  write(std_out,*)' critic : exit, ires=1'
319 !  ENDDEBUG
320    return
321 !  stop 'ERROR INVERTING HESSIAN'
322  end if
323  do ii=1,3
324    yy(ii,1:3)=0.
325    yy(ii,ii)=1.
326  end do
327  do jii=1,3
328    call lubksb(hrho,3,3,ipiv,yy(1,jii))
329  end do
330 
331 
332 !write(std_out,'(":HESSIAN^(-1) ",/,3F16.8,/,3F16.8,/,3F16.8)') ((y(ii,jii),jii=1,3),ii=1,3)
333 
334  call vgh_rho(vv,chg,grho,hrho,rr,iat,ipos,0)
335 
336 !write(std_out,'("LAPLAC:",F16.8)') hrho(1,1)+hrho(2,2)+hrho(3,3)
337 
338  call jacobi(hrho,3,3,ev,yy,nrot)
339  call ordr(ev,yy,3,1)
340  zz(:,:)=yy(:,:)
341 
342 !do ii=1,3
343 !do jii=1,3
344 !zz(ii,jii)=yy(jii,ii)
345 !end do
346 !end do
347 
348 !write(std_out,'(":AUTOVAL ",3F16.8)') (ev(ii),ii=1,3)
349 !write(std_out,'(":AUTOVEC ",/,3F16.8,/,3F16.8,/,3F16.8)') ((zz(ii,jii),ii=1,3),jii=1,3)
350 
351  if (sort/=0)  then
352    ABI_DEALLOCATE(pom)
353    ABI_DEALLOCATE(pom2)
354    ABI_DEALLOCATE(lamb)
355  end if
356 
357 !DEBUG
358 !write(std_out,*)' critic : exit, ires= ',ires
359 !ENDDEBUG
360 end subroutine critic

ABINIT/ordr [ Functions ]

[ Top ] [ Functions ]

NAME

 ordr

FUNCTION

COPYRIGHT

 Copyright (C) 2007-2017 ABINIT group ( ).
 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

  (to be filled)

OUTPUT

  (to be filled)

PARENTS

      critic

CHILDREN

SOURCE

391 #if defined HAVE_CONFIG_H
392 #include "config.h"
393 #endif
394 
395 subroutine ordr(aa,dd,nn,cff)
396 
397  use m_profiling_abi
398 
399  use defs_basis
400 
401 !This section has been created automatically by the script Abilint (TD).
402 !Do not modify the following lines by hand.
403 #undef ABI_FUNC
404 #define ABI_FUNC 'ordr'
405 !End of the abilint section
406 
407  implicit none
408 
409 !Arguments ----------------------------
410 !scalars
411  integer,intent(in) :: cff,nn
412 !arrays
413  real(dp),intent(inout) :: aa(nn),dd(nn,nn)
414 
415 !Local variables ----------------------
416 !scalars
417  integer :: ii,jj,kk
418  real(dp) :: uu
419 
420 ! *********************************************************************
421 
422  do ii=1,nn-1
423    kk=ii
424    uu=aa(ii)
425    do jj=ii+1,nn
426      if (cff==1) then
427        if (aa(jj) >= uu+tol12) then
428          kk=jj
429          uu=aa(jj)
430        end if
431      else
432        if (aa(jj) <= uu-tol12) then
433          kk=jj
434          uu=aa(jj)
435        end if
436      end if
437    end do
438    if (kk /= ii) then
439      aa(kk)=aa(ii)
440      aa(ii)=uu
441      do jj=1,nn
442        uu=dd(jj,ii)
443        dd(jj,ii)=dd(jj,kk)
444        dd(jj,kk)=uu
445      end do
446    end if
447  end do
448 end subroutine ordr