TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
NOTES
PARENTS
mover
CHILDREN
bonds_free,dsyev,dsysv,fcart2fred,make_bonds_new,xred2xcart
SOURCE
40 #if defined HAVE_CONFIG_H 41 #include "config.h" 42 #endif 43 44 #include "abi_common.h" 45 46 subroutine prec_simple(ab_mover,forstr,hist,icycle,itime,iexit) 47 48 use defs_basis 49 use m_profiling_abi 50 use m_abimover 51 use m_abihist 52 use m_linalg_interfaces 53 54 use m_geometry, only : fcart2fred, xred2xcart 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 'prec_simple' 60 !End of the abilint section 61 62 implicit none 63 64 !Arguments ------------------------------------ 65 !scalars 66 integer,intent(in) :: iexit,itime,icycle 67 type(abimover),intent(in) :: ab_mover 68 type(abihist),intent(inout) :: hist 69 type(abiforstr),intent(inout) :: forstr 70 71 !Local variables------------------------------- 72 !scalars 73 integer :: period,ii,jj,index,kk,ksub,jsub 74 integer :: info,lwork,new_order_forces 75 real(dp) :: Z,badgerfactor,lambda,sigma,val_rms 76 integer,save :: order_forces 77 logical :: Compute_Matrix 78 !arrays 79 type(go_bonds) :: bonds 80 integer,allocatable :: periods(:,:) 81 integer,allocatable :: iatoms(:,:) 82 integer :: ipiv(3*ab_mover%natom) 83 real(dp) :: xcart(3,ab_mover%natom) 84 real(dp) :: fcart(3,ab_mover%natom) 85 real(dp) :: B(3*ab_mover%natom) 86 real(dp) :: rprimd(3,3) 87 real(dp) :: matrix_tmp(3*ab_mover%natom,3*ab_mover%natom) 88 real(dp) :: w(3*ab_mover%natom) 89 real(dp),allocatable :: work(:) 90 real(dp) :: badger(6,6) 91 real(dp),allocatable,save :: matrix(:,:) 92 character(len=18) :: fmt 93 94 !*************************************************************************** 95 !Beginning of executable session 96 !*************************************************************************** 97 98 if (iexit/=0)then 99 ABI_DEALLOCATE(matrix) 100 return 101 end if 102 103 !########################################################## 104 !### 01. Show the Precondition parameters, set the badger 105 !### matrix. 106 107 write(std_out,*) 'Precondition option',ab_mover%goprecon 108 write(std_out,*) 'Precondition parameters',ab_mover%goprecprm 109 lambda=ab_mover%goprecprm(1) 110 111 badger(:,:)=reshape( (/ -0.2573, 0.3401, 0.6937, 0.7126, 0.8335, 0.9491,& 112 & 0.3401, 0.9652, 1.2843, 1.4725, 1.6549, 1.7190,& 113 & 0.6937, 1.2843, 1.6925, 1.8238, 2.1164, 2.3185,& 114 & 0.7126, 1.4725, 1.8238, 2.0203, 2.2137, 2.5206,& 115 & 0.8335, 1.6549, 2.1164, 2.2137, 2.3718, 2.5110,& 116 & 0.9491, 1.7190, 2.3185, 2.5206, 2.5110, 0.0000 /), (/ 6, 6/) ) 117 118 write(fmt,'(a1,i4,a5)') '(',3*ab_mover%natom,'f8.3)' 119 120 !########################################################## 121 !### 02. Take the coordinates and cell parameters from HIST 122 123 rprimd(:,:)=hist%rprimd(:,:,hist%ihist) 124 fcart(:,:)=hist%fcart(:,:,hist%ihist) 125 call xred2xcart(ab_mover%natom,rprimd,xcart,hist%xred(:,:,hist%ihist)) 126 127 !########################################################## 128 !### 03. Decide based on kind of precondiotioner if 129 !### a new matrix should be computed 130 131 if (ab_mover%goprecon==2)then 132 133 val_rms=0.0 134 do kk=1,ab_mover%natom 135 do jj=1,3 136 val_rms=val_rms+fcart(jj,kk)**2 137 end do 138 end do 139 val_rms=sqrt(val_rms/dble(ab_mover%natom)) 140 new_order_forces=int(log(val_rms)/log(10.0)) 141 end if 142 143 if (itime==1.and.icycle==1)then 144 Compute_Matrix=.TRUE. 145 order_forces=new_order_forces 146 if (allocated(matrix)) then 147 ABI_DEALLOCATE(matrix) 148 end if 149 150 ABI_ALLOCATE(matrix,(3*ab_mover%natom,3*ab_mover%natom)) 151 else 152 Compute_Matrix=.FALSE. 153 if ((ab_mover%goprecon==2).and.(order_forces.gt.new_order_forces)) then 154 Compute_Matrix=.TRUE. 155 order_forces=new_order_forces 156 end if 157 if (ab_mover%goprecon==3) Compute_Matrix=.TRUE. 158 end if 159 160 !########################################################## 161 !### 04. Compute a new precondition matrix if required 162 163 if (Compute_Matrix)then 164 165 ! Fix the tolerance for create a bond 166 bonds%tolerance=1.35 167 bonds%nbonds=1 168 169 ! Allocate the arrays with exactly the rigth nbonds 170 ABI_ALLOCATE(bonds%bond_vect,(3,bonds%nbonds)) 171 ABI_ALLOCATE(bonds%bond_length,(bonds%nbonds)) 172 ABI_ALLOCATE(bonds%indexi,(ab_mover%natom,bonds%nbonds)) 173 ABI_ALLOCATE(bonds%nbondi,(ab_mover%natom)) 174 175 ! Compute the bonds 176 call make_bonds_new(bonds,ab_mover%natom,ab_mover%ntypat,rprimd,& 177 & ab_mover%typat,xcart,ab_mover%znucl) 178 179 write(std_out,'(a,a,63a,a)') ch10,'---PRECONDITIONER',('-',kk=1,63),ch10 180 181 ! For all bonds detect wich atoms are involved 182 ! and wich period they coprrespond in the periodic table 183 if (bonds%nbonds>0)then 184 185 ABI_ALLOCATE(periods,(2,bonds%nbonds)) 186 ABI_ALLOCATE(iatoms,(2,bonds%nbonds)) 187 periods(:,:)=0 188 iatoms(:,:)=0 189 190 write(std_out,'(a)') 'Bond of Atom | Bond Number | Index' 191 192 do ii=1,ab_mover%natom 193 Z=ab_mover%znucl(ab_mover%typat(ii)) 194 if (Z==1)then 195 period=1 196 elseif ((Z>1).and.(Z<10))then 197 period=2 198 elseif ((Z>10).and.(Z<18))then 199 period=3 200 elseif ((Z>18).and.(Z<36))then 201 period=4 202 elseif ((Z>36).and.(Z<54))then 203 period=5 204 elseif ((Z>55).and.(Z<86))then 205 period=6 206 else 207 ! Here are the cases for atoms larger than Fr(87) and 208 ! All the noble gases He-Rn 209 period=-1 210 end if 211 212 do jj=1,bonds%nbondi(ii) 213 index=bonds%indexi(ii,jj) 214 215 write(std_out,'(i6,a,i6,a,i4)') ii,' |',jj,' |',index 216 217 ! The first atom should have index=0 218 ! To make easy fill the matrix using its 219 ! index 220 221 if (index>0)then 222 periods(1,index)=period 223 iatoms(1,index)=ii 224 elseif (index<0) then 225 periods(2,-index)=period 226 iatoms(2,-index)=ii 227 end if 228 end do 229 end do 230 231 write(std_out,'(a)') ch10 232 233 end if 234 235 ! For all bonds compute the 3x3 matrix and fill also the big matrix 236 237 matrix(:,:)=0.0_dp 238 do ii=1,bonds%nbonds 239 240 write(std_out,*) 'Bond number:',ii 241 if (iatoms(1,ii)>0 .and. iatoms(2,ii)>0) then 242 write(std_out,*) 'Between atoms:',iatoms(1,ii),' and ',iatoms(2,ii) 243 badgerfactor=badger(periods(1,ii),periods(2,ii)) 244 write(std_out,*) 'Periods of atoms:',periods(1,ii),' and ',periods(2,ii) 245 write(std_out,*) 'Badger factor:',badgerfactor 246 247 ! Compute the diadic product and 248 ! Insert the matrix into the big one 249 do jj=1,3 250 do kk=1,3 251 ! The non diagonal elements 252 jsub=3*(iatoms(1,ii)-1)+jj 253 ksub=3*(iatoms(2,ii)-1)+kk 254 matrix(jsub,ksub)=matrix(jsub,ksub)-& 255 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 256 257 jsub=3*(iatoms(2,ii)-1)+jj 258 ksub=3*(iatoms(1,ii)-1)+kk 259 matrix(jsub,ksub)=matrix(jsub,ksub)-& 260 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 261 262 ! The diagonal blocks 263 jsub=3*(iatoms(1,ii)-1)+jj 264 ksub=3*(iatoms(1,ii)-1)+kk 265 matrix(jsub,ksub)=matrix(jsub,ksub)+& 266 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 267 268 jsub=3*(iatoms(2,ii)-1)+jj 269 ksub=3*(iatoms(2,ii)-1)+kk 270 matrix(jsub,ksub)=matrix(jsub,ksub)+& 271 & badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii) 272 273 end do !do kk=1,3 274 end do !do jj=1,3 275 276 end if 277 278 end do 279 280 if (bonds%nbonds>0)then 281 ABI_DEALLOCATE(periods) 282 ABI_DEALLOCATE(iatoms) 283 end if 284 285 call bonds_free(bonds) 286 287 if (3*ab_mover%natom<100)then 288 ! Visualize the matrix 289 do jj=1,3*ab_mover%natom 290 write (std_out,fmt) matrix(jj,:) 291 end do 292 end if 293 294 matrix_tmp(:,:)=matrix(:,:) 295 !write(*,*)"matrix_tmp",matrix_tmp 296 297 ABI_ALLOCATE(work,(1)) 298 lwork=-1 299 call DSYEV('V', 'U', 3*ab_mover%natom, matrix_tmp, 3*ab_mover%natom, w , work, lwork, info ) 300 lwork=work(1) 301 write(std_out,*) '[DSYEV] Recommended lwork=',lwork 302 ABI_DEALLOCATE(work) 303 ABI_ALLOCATE(work,(lwork)) 304 call DSYEV('V', 'U', 3*ab_mover%natom, matrix_tmp, 3*ab_mover%natom, w , work, lwork, info ) 305 ABI_DEALLOCATE(work) 306 307 write(std_out,*) 'DSYEV info=',info 308 write(std_out,*) 'Eigenvalues:' 309 write(std_out,fmt) w(:) 310 311 sigma=0 312 do jj=1,3*ab_mover%natom 313 sigma=max(w(jj),sigma) 314 end do 315 316 matrix=lambda*matrix 317 318 write(std_out,*) ch10 319 do ii=1,3*ab_mover%natom 320 matrix(ii,ii)=matrix(ii,ii)+(1-lambda)*sigma 321 end do 322 323 end if ! if (Compute_Matrix) 324 325 !########################################################## 326 !### 05. Use the precondition matrix to compute new residuals 327 328 B=reshape(fcart,(/ 3*ab_mover%natom /)) 329 330 if (3*ab_mover%natom<100)then 331 ! Visualize the matrix 332 do jj=1,3*ab_mover%natom 333 write (std_out,fmt) matrix(jj,:) 334 end do 335 end if 336 337 !call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info ) 338 !MGNAG FIXME: This call causes a floating point exception if NAG+MKL 339 ABI_ALLOCATE(work,(1)) 340 lwork=-1 341 call DSYSV( 'U', 3*ab_mover%natom, 1, matrix,& 342 & 3*ab_mover%natom, ipiv, B, 3*ab_mover%natom, work, lwork, info ) 343 344 lwork=work(1) 345 write(std_out,*) '[DSYSV] Recomended lwork=',lwork 346 ABI_DEALLOCATE(work) 347 ABI_ALLOCATE(work,(lwork)) 348 call DSYSV( 'U', 3*ab_mover%natom, 1, matrix,& 349 & 3*ab_mover%natom, ipiv, B, 3*ab_mover%natom, work, lwork, info ) 350 ABI_DEALLOCATE(work) 351 352 write(std_out,*) 'DSYSV info=',info 353 write(std_out,*) 'Solution:' 354 write(std_out,fmt) B(:) 355 356 forstr%fcart=reshape(B,(/ 3, ab_mover%natom /) ) 357 call fcart2fred(forstr%fcart,forstr%fred,rprimd,ab_mover%natom) 358 359 end subroutine prec_simple