TABLE OF CONTENTS
ABINIT/critic [ 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 ]
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