TABLE OF CONTENTS
ABINIT/m_pred_simple [ Modules ]
NAME
m_pred_simple
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, JCC, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_pred_simple 23 24 use defs_basis 25 use m_abicore 26 use m_abimover 27 use m_abihist 28 29 use m_geometry, only : fcart2gred, xred2xcart 30 31 implicit none 32 33 private
ABINIT/prec_simple [ Functions ]
NAME
prec_simple
FUNCTION
Simple preconditioner, compute the force constant matrix
using the Badger's rule:
F=A/(r-B)^3
INPUTS
ab_mover <type(abimover)> : Datatype with all the information needed by the preconditioner zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
142 subroutine prec_simple(ab_mover,forstr,hist,icycle,itime,iexit) 143 144 use m_linalg_interfaces 145 146 !Arguments ------------------------------------ 147 !scalars 148 integer,intent(in) :: iexit,itime,icycle 149 type(abimover),intent(in) :: ab_mover 150 type(abihist),intent(in) :: hist 151 type(abiforstr),intent(inout) :: forstr 152 153 !Local variables------------------------------- 154 !scalars 155 integer :: period,ii,jj,index,kk,ksub,jsub 156 integer :: info,lwork,new_order_forces 157 real(dp) :: Z,badgerfactor,lambda,sigma,val_rms 158 integer,save :: order_forces 159 logical :: Compute_Matrix 160 !arrays 161 type(go_bonds) :: bonds 162 integer,allocatable :: periods(:,:) 163 integer,allocatable :: iatoms(:,:) 164 integer :: ipiv(3*ab_mover%natom) 165 real(dp) :: xcart(3,ab_mover%natom) 166 real(dp) :: fcart(3,ab_mover%natom) 167 real(dp) :: B(3*ab_mover%natom) 168 real(dp) :: rprimd(3,3) 169 real(dp) :: w(3*ab_mover%natom) 170 real(dp),allocatable :: matrix_tmp(:,:) 171 real(dp),allocatable :: work(:) 172 real(dp) :: badger(6,6) 173 real(dp),allocatable,save :: matrix(:,:) 174 character(len=18) :: fmt 175 176 !*************************************************************************** 177 !Beginning of executable session 178 !*************************************************************************** 179 180 if (iexit/=0)then 181 if(allocated(matrix))then 182 ABI_FREE(matrix) 183 endif 184 return 185 end if 186 187 !########################################################## 188 !### 01. Show the Precondition parameters, set the badger 189 !### matrix. 190 191 write(std_out,*) 'Precondition option',ab_mover%goprecon 192 write(std_out,*) 'Precondition parameters',ab_mover%goprecprm 193 lambda=ab_mover%goprecprm(1) 194 195 badger(:,:)=reshape( (/ -0.2573, 0.3401, 0.6937, 0.7126, 0.8335, 0.9491,& 196 & 0.3401, 0.9652, 1.2843, 1.4725, 1.6549, 1.7190,& 197 & 0.6937, 1.2843, 1.6925, 1.8238, 2.1164, 2.3185,& 198 & 0.7126, 1.4725, 1.8238, 2.0203, 2.2137, 2.5206,& 199 & 0.8335, 1.6549, 2.1164, 2.2137, 2.3718, 2.5110,& 200 & 0.9491, 1.7190, 2.3185, 2.5206, 2.5110, 0.0000 /), (/ 6, 6/) ) 201 202 write(fmt,'(a1,i4,a5)') '(',3*ab_mover%natom,'f8.3)' 203 204 !########################################################## 205 !### 02. Take the coordinates and cell parameters from HIST 206 207 rprimd(:,:)=hist%rprimd(:,:,hist%ihist) 208 fcart(:,:)=hist%fcart(:,:,hist%ihist) 209 call xred2xcart(ab_mover%natom,rprimd,xcart,hist%xred(:,:,hist%ihist)) 210 211 !########################################################## 212 !### 03. Decide based on kind of preconditioner if 213 !### a new matrix should be computed 214 215 new_order_forces = one ! This to avoid using unitialized variables. 216 217 if (ab_mover%goprecon==2)then 218 219 val_rms=0.0 220 do kk=1,ab_mover%natom 221 do jj=1,3 222 val_rms=val_rms+fcart(jj,kk)**2 223 end do 224 end do 225 val_rms=sqrt(val_rms/dble(ab_mover%natom)) 226 new_order_forces=int(log(val_rms)/log(10.0)) 227 end if 228 229 if (itime==1.and.icycle==1)then 230 Compute_Matrix=.TRUE. 231 order_forces=new_order_forces 232 if (allocated(matrix)) then 233 ABI_FREE(matrix) 234 end if 235 236 ABI_MALLOC(matrix,(3*ab_mover%natom,3*ab_mover%natom)) 237 else 238 Compute_Matrix=.FALSE. 239 if ((ab_mover%goprecon==2).and.(order_forces.gt.new_order_forces)) then 240 Compute_Matrix=.TRUE. 241 order_forces=new_order_forces 242 end if 243 if (ab_mover%goprecon==3) Compute_Matrix=.TRUE. 244 end if 245 246 !########################################################## 247 !### 04. Compute a new precondition matrix if required 248 249 if (Compute_Matrix)then 250 251 ! Fix the tolerance for create a bond 252 bonds%tolerance=1.35 253 bonds%nbonds=1 254 255 ! Allocate the arrays with exactly the rigth nbonds 256 ABI_MALLOC(bonds%bond_vect,(3,bonds%nbonds)) 257 ABI_MALLOC(bonds%bond_length,(bonds%nbonds)) 258 ABI_MALLOC(bonds%indexi,(ab_mover%natom,bonds%nbonds)) 259 ABI_MALLOC(bonds%nbondi,(ab_mover%natom)) 260 261 ! Compute the bonds 262 call make_bonds_new(bonds,ab_mover%natom,ab_mover%ntypat,rprimd,& 263 & ab_mover%typat,xcart,ab_mover%znucl) 264 265 write(std_out,'(a,a,63a,a)') ch10,'---PRECONDITIONER',('-',kk=1,63),ch10 266 267 ! For all bonds detect wich atoms are involved 268 ! and wich period they coprrespond in the periodic table 269 if (bonds%nbonds>0)then 270 271 ABI_MALLOC(periods,(2,bonds%nbonds)) 272 ABI_MALLOC(iatoms,(2,bonds%nbonds)) 273 periods(:,:)=0 274 iatoms(:,:)=0 275 276 write(std_out,'(a)') 'Bond of Atom | Bond Number | Index' 277 278 do ii=1,ab_mover%natom 279 Z=ab_mover%znucl(ab_mover%typat(ii)) 280 if (Z==1)then 281 period=1 282 elseif ((Z>1).and.(Z<10))then 283 period=2 284 elseif ((Z>10).and.(Z<18))then 285 period=3 286 elseif ((Z>18).and.(Z<36))then 287 period=4 288 elseif ((Z>36).and.(Z<54))then 289 period=5 290 elseif ((Z>55).and.(Z<86))then 291 period=6 292 else 293 ! Here are the cases for atoms larger than Fr(87) and 294 ! All the noble gases He-Rn 295 period=-1 296 end if 297 298 do jj=1,bonds%nbondi(ii) 299 index=bonds%indexi(ii,jj) 300 301 write(std_out,'(i6,a,i6,a,i4)') ii,' |',jj,' |',index 302 303 ! The first atom should have index=0 304 ! To make easy fill the matrix using its 305 ! index 306 307 if (index>0)then 308 periods(1,index)=period 309 iatoms(1,index)=ii 310 elseif (index<0) then 311 periods(2,-index)=period 312 iatoms(2,-index)=ii 313 end if 314 end do 315 end do 316 317 write(std_out,'(a)') ch10 318 319 end if 320 321 ! For all bonds compute the 3x3 matrix and fill also the big matrix 322 323 matrix(:,:)=0.0_dp 324 do ii=1,bonds%nbonds 325 326 write(std_out,*) 'Bond number:',ii 327 if (iatoms(1,ii)>0 .and. iatoms(2,ii)>0) then 328 write(std_out,*) 'Between atoms:',iatoms(1,ii),' and ',iatoms(2,ii) 329 badgerfactor=badger(periods(1,ii),periods(2,ii)) 330 write(std_out,*) 'Periods of atoms:',periods(1,ii),' and ',periods(2,ii) 331 write(std_out,*) 'Badger factor:',badgerfactor 332 333 ! Compute the diadic product and 334 ! Insert the matrix into the big one 335 do jj=1,3 336 do kk=1,3 337 ! The non diagonal elements 338 jsub=3*(iatoms(1,ii)-1)+jj 339 ksub=3*(iatoms(2,ii)-1)+kk 340 matrix(jsub,ksub)=matrix(jsub,ksub)-& 341 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 342 343 jsub=3*(iatoms(2,ii)-1)+jj 344 ksub=3*(iatoms(1,ii)-1)+kk 345 matrix(jsub,ksub)=matrix(jsub,ksub)-& 346 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 347 348 ! The diagonal blocks 349 jsub=3*(iatoms(1,ii)-1)+jj 350 ksub=3*(iatoms(1,ii)-1)+kk 351 matrix(jsub,ksub)=matrix(jsub,ksub)+& 352 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 353 354 jsub=3*(iatoms(2,ii)-1)+jj 355 ksub=3*(iatoms(2,ii)-1)+kk 356 matrix(jsub,ksub)=matrix(jsub,ksub)+& 357 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 358 359 end do !do kk=1,3 360 end do !do jj=1,3 361 362 end if 363 364 end do 365 366 if (bonds%nbonds>0)then 367 ABI_FREE(periods) 368 ABI_FREE(iatoms) 369 end if 370 371 call bonds_free(bonds) 372 373 if (3*ab_mover%natom<100)then 374 ! Visualize the matrix 375 do jj=1,3*ab_mover%natom 376 write (std_out,fmt) matrix(jj,:) 377 end do 378 end if 379 380 ABI_MALLOC(matrix_tmp,(3*ab_mover%natom,3*ab_mover%natom)) 381 382 matrix_tmp(:,:)=matrix(:,:) 383 !write(*,*)"matrix_tmp",matrix_tmp 384 385 ABI_MALLOC(work,(1)) 386 lwork=-1 387 call DSYEV('V', 'U', 3*ab_mover%natom, matrix_tmp, 3*ab_mover%natom, w , work, lwork, info ) 388 lwork=work(1) 389 write(std_out,*) '[DSYEV] Recommended lwork=',lwork 390 ABI_FREE(work) 391 ABI_MALLOC(work,(lwork)) 392 call DSYEV('V', 'U', 3*ab_mover%natom, matrix_tmp, 3*ab_mover%natom, w , work, lwork, info ) 393 ABI_FREE(work) 394 ABI_FREE(matrix_tmp) 395 396 write(std_out,*) 'DSYEV info=',info 397 write(std_out,*) 'Eigenvalues:' 398 write(std_out,fmt) w(:) 399 400 sigma=0 401 do jj=1,3*ab_mover%natom 402 sigma=max(w(jj),sigma) 403 end do 404 405 matrix=lambda*matrix 406 407 write(std_out,*) ch10 408 do ii=1,3*ab_mover%natom 409 matrix(ii,ii)=matrix(ii,ii)+(1-lambda)*sigma 410 end do 411 412 end if ! if (Compute_Matrix) 413 414 !########################################################## 415 !### 05. Use the precondition matrix to compute new residuals 416 417 B=reshape(fcart,(/ 3*ab_mover%natom /)) 418 419 if (3*ab_mover%natom<100)then 420 ! Visualize the matrix 421 do jj=1,3*ab_mover%natom 422 write (std_out,fmt) matrix(jj,:) 423 end do 424 end if 425 426 !call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info ) 427 !MGNAG FIXME: This call causes a floating point exception if NAG+MKL 428 ABI_MALLOC(work,(1)) 429 lwork=-1 430 call DSYSV( 'U', 3*ab_mover%natom, 1, matrix,& 431 & 3*ab_mover%natom, ipiv, B, 3*ab_mover%natom, work, lwork, info ) 432 433 lwork=work(1) 434 write(std_out,*) '[DSYSV] Recomended lwork=',lwork 435 ABI_FREE(work) 436 ABI_MALLOC(work,(lwork)) 437 call DSYSV( 'U', 3*ab_mover%natom, 1, matrix,& 438 & 3*ab_mover%natom, ipiv, B, 3*ab_mover%natom, work, lwork, info ) 439 ABI_FREE(work) 440 441 write(std_out,*) 'DSYSV info=',info 442 write(std_out,*) 'Solution:' 443 write(std_out,fmt) B(:) 444 445 forstr%fcart=reshape(B,(/ 3, ab_mover%natom /) ) 446 call fcart2gred(forstr%fcart,forstr%gred,rprimd,ab_mover%natom) 447 448 end subroutine prec_simple
ABINIT/pred_simple [ Functions ]
NAME
pred_simple
FUNCTION
Ionmov predictors (4 & 5) Internal to scfcv. Actually, this routine does nothing (only copy) as all operations are internal to scfcv ... IONMOV 4: Conjugate gradient algorithm for simultaneous optimization of potential and ionic degrees of freedom. It can be used with iscf=2 and iscf=5 or 6 IONMOV 5: Simple relaxation of ionic positions according to (converged) forces. Equivalent to ionmov=1 with zero masses, albeit the relaxation coefficient is not vis, but iprcfc.
INPUTS
ab_mover <type(abimover)> : Datatype with all the information needed by the preditor zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
72 subroutine pred_simple(ab_mover,hist,iexit) 73 74 !Arguments ------------------------------------ 75 !scalars 76 type(abimover),intent(in) :: ab_mover 77 type(abihist),intent(inout) :: hist 78 integer,intent(in) :: iexit 79 80 !Local variables------------------------------- 81 !scalars 82 integer :: ihist_next,kk,jj 83 84 !*************************************************************************** 85 !Beginning of executable session 86 !*************************************************************************** 87 88 if(iexit/=0)then 89 return 90 end if 91 92 !All the operations are internal to scfcv.F90 93 94 !XRED, FCART and VEL 95 ihist_next = abihist_findIndex(hist,+1) 96 do kk=1,ab_mover%natom 97 do jj=1,3 98 hist%xred(jj,kk, ihist_next)=hist%xred (jj,kk,hist%ihist) 99 hist%fcart(jj,kk,ihist_next)=hist%fcart(jj,kk,hist%ihist) 100 hist%vel(jj,kk,ihist_next)=hist%vel(jj,kk,hist%ihist) 101 end do 102 end do 103 104 !ACELL 105 do jj=1,3 106 hist%acell(jj,ihist_next)=hist%acell(jj,hist%ihist) 107 end do 108 109 !RPRIMD 110 do kk=1,3 111 do jj=1,3 112 hist%rprimd(jj,kk,ihist_next)=hist%rprimd(jj,kk,hist%ihist) 113 end do 114 end do 115 116 hist%ihist=ihist_next 117 118 end subroutine pred_simple