TABLE OF CONTENTS


ABINIT/m_pred_simple [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_simple

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_pred_simple
28 
29  use defs_basis
30  use m_abicore
31  use m_abimover
32  use m_abihist
33 
34  use m_geometry,  only : fcart2fred, xred2xcart
35 
36  implicit none
37 
38  private

ABINIT/prec_simple [ Functions ]

[ Top ] [ 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

PARENTS

      mover

CHILDREN

      bonds_free,dsyev,dsysv,fcart2fred,make_bonds_new,xred2xcart

SOURCE

166 subroutine prec_simple(ab_mover,forstr,hist,icycle,itime,iexit)
167 
168  use m_linalg_interfaces
169 
170 !This section has been created automatically by the script Abilint (TD).
171 !Do not modify the following lines by hand.
172 #undef ABI_FUNC
173 #define ABI_FUNC 'prec_simple'
174 !End of the abilint section
175 
176  implicit none
177 
178 !Arguments ------------------------------------
179 !scalars
180  integer,intent(in) :: iexit,itime,icycle
181  type(abimover),intent(in) :: ab_mover
182  type(abihist),intent(in) :: hist
183  type(abiforstr),intent(inout) :: forstr
184 
185 !Local variables-------------------------------
186 !scalars
187  integer :: period,ii,jj,index,kk,ksub,jsub
188  integer :: info,lwork,new_order_forces
189  real(dp) :: Z,badgerfactor,lambda,sigma,val_rms
190  integer,save :: order_forces
191  logical :: Compute_Matrix
192 !arrays
193  type(go_bonds) :: bonds
194  integer,allocatable :: periods(:,:)
195  integer,allocatable :: iatoms(:,:)
196  integer  :: ipiv(3*ab_mover%natom)
197  real(dp) :: xcart(3,ab_mover%natom)
198  real(dp) :: fcart(3,ab_mover%natom)
199  real(dp) :: B(3*ab_mover%natom)
200  real(dp) :: rprimd(3,3)
201  real(dp) :: w(3*ab_mover%natom)
202  real(dp),allocatable :: matrix_tmp(:,:)
203  real(dp),allocatable :: work(:)
204  real(dp) :: badger(6,6)
205  real(dp),allocatable,save :: matrix(:,:)
206  character(len=18)   :: fmt
207 
208 !***************************************************************************
209 !Beginning of executable session
210 !***************************************************************************
211 
212  if (iexit/=0)then
213    if(allocated(matrix))then
214      ABI_DEALLOCATE(matrix)
215    endif
216    return
217  end if
218 
219 !##########################################################
220 !### 01. Show the Precondition parameters, set the badger
221 !###     matrix.
222 
223  write(std_out,*) 'Precondition option',ab_mover%goprecon
224  write(std_out,*) 'Precondition parameters',ab_mover%goprecprm
225  lambda=ab_mover%goprecprm(1)
226 
227  badger(:,:)=reshape( (/ -0.2573, 0.3401, 0.6937, 0.7126, 0.8335, 0.9491,&
228 & 0.3401, 0.9652, 1.2843, 1.4725, 1.6549, 1.7190,&
229 & 0.6937, 1.2843, 1.6925, 1.8238, 2.1164, 2.3185,&
230 & 0.7126, 1.4725, 1.8238, 2.0203, 2.2137, 2.5206,&
231 & 0.8335, 1.6549, 2.1164, 2.2137, 2.3718, 2.5110,&
232 & 0.9491, 1.7190, 2.3185, 2.5206, 2.5110, 0.0000 /), (/ 6, 6/) )
233 
234  write(fmt,'(a1,i4,a5)') '(',3*ab_mover%natom,'f8.3)'
235 
236 !##########################################################
237 !### 02. Take the coordinates and cell parameters from HIST
238 
239  rprimd(:,:)=hist%rprimd(:,:,hist%ihist)
240  fcart(:,:)=hist%fcart(:,:,hist%ihist)
241  call xred2xcart(ab_mover%natom,rprimd,xcart,hist%xred(:,:,hist%ihist))
242 
243 !##########################################################
244 !### 03. Decide based on kind of preconditioner if
245 !###     a new matrix should be computed
246 
247  if (ab_mover%goprecon==2)then
248 
249    val_rms=0.0
250    do kk=1,ab_mover%natom
251      do jj=1,3
252        val_rms=val_rms+fcart(jj,kk)**2
253      end do
254    end do
255    val_rms=sqrt(val_rms/dble(ab_mover%natom))
256    new_order_forces=int(log(val_rms)/log(10.0))
257  end if
258 
259  if (itime==1.and.icycle==1)then
260    Compute_Matrix=.TRUE.
261    order_forces=new_order_forces
262    if (allocated(matrix))  then
263      ABI_DEALLOCATE(matrix)
264    end if
265 
266    ABI_ALLOCATE(matrix,(3*ab_mover%natom,3*ab_mover%natom))
267  else
268    Compute_Matrix=.FALSE.
269    if ((ab_mover%goprecon==2).and.(order_forces.gt.new_order_forces)) then
270      Compute_Matrix=.TRUE.
271      order_forces=new_order_forces
272    end if
273    if (ab_mover%goprecon==3) Compute_Matrix=.TRUE.
274  end if
275 
276 !##########################################################
277 !### 04. Compute a new precondition matrix if required
278 
279  if (Compute_Matrix)then
280 
281 !  Fix the tolerance for create a bond
282    bonds%tolerance=1.35
283    bonds%nbonds=1
284 
285 !  Allocate the arrays with exactly the rigth nbonds
286    ABI_ALLOCATE(bonds%bond_vect,(3,bonds%nbonds))
287    ABI_ALLOCATE(bonds%bond_length,(bonds%nbonds))
288    ABI_ALLOCATE(bonds%indexi,(ab_mover%natom,bonds%nbonds))
289    ABI_ALLOCATE(bonds%nbondi,(ab_mover%natom))
290 
291 !  Compute the bonds
292    call make_bonds_new(bonds,ab_mover%natom,ab_mover%ntypat,rprimd,&
293 &   ab_mover%typat,xcart,ab_mover%znucl)
294 
295    write(std_out,'(a,a,63a,a)') ch10,'---PRECONDITIONER',('-',kk=1,63),ch10
296 
297 !  For all bonds detect wich atoms are involved
298 !  and wich period they coprrespond in the periodic table
299    if (bonds%nbonds>0)then
300 
301      ABI_ALLOCATE(periods,(2,bonds%nbonds))
302      ABI_ALLOCATE(iatoms,(2,bonds%nbonds))
303      periods(:,:)=0
304      iatoms(:,:)=0
305 
306      write(std_out,'(a)') 'Bond of Atom | Bond Number | Index'
307 
308      do ii=1,ab_mover%natom
309        Z=ab_mover%znucl(ab_mover%typat(ii))
310        if (Z==1)then
311          period=1
312        elseif ((Z>1).and.(Z<10))then
313          period=2
314        elseif ((Z>10).and.(Z<18))then
315          period=3
316        elseif ((Z>18).and.(Z<36))then
317          period=4
318        elseif ((Z>36).and.(Z<54))then
319          period=5
320        elseif ((Z>55).and.(Z<86))then
321          period=6
322        else
323 !        Here are the cases for atoms larger than Fr(87) and
324 !        All the noble gases He-Rn
325          period=-1
326        end if
327 
328        do jj=1,bonds%nbondi(ii)
329          index=bonds%indexi(ii,jj)
330 
331          write(std_out,'(i6,a,i6,a,i4)') ii,'       |',jj,'       |',index
332 
333 !        The first atom should have index=0
334 !        To make easy fill the matrix using its
335 !        index
336 
337          if (index>0)then
338            periods(1,index)=period
339            iatoms(1,index)=ii
340          elseif (index<0) then
341            periods(2,-index)=period
342            iatoms(2,-index)=ii
343          end if
344        end do
345      end do
346 
347      write(std_out,'(a)') ch10
348 
349    end if
350 
351 !  For all bonds compute the 3x3 matrix and fill also the big matrix
352 
353    matrix(:,:)=0.0_dp
354    do ii=1,bonds%nbonds
355 
356      write(std_out,*) 'Bond number:',ii
357      if (iatoms(1,ii)>0 .and. iatoms(2,ii)>0) then
358        write(std_out,*) 'Between atoms:',iatoms(1,ii),' and ',iatoms(2,ii)
359        badgerfactor=badger(periods(1,ii),periods(2,ii))
360        write(std_out,*) 'Periods of atoms:',periods(1,ii),' and ',periods(2,ii)
361        write(std_out,*) 'Badger factor:',badgerfactor
362 
363 !      Compute the diadic product and
364 !      Insert the matrix into the big one
365        do jj=1,3
366          do kk=1,3
367 !          The non diagonal elements
368            jsub=3*(iatoms(1,ii)-1)+jj
369            ksub=3*(iatoms(2,ii)-1)+kk
370            matrix(jsub,ksub)=matrix(jsub,ksub)-&
371 &           badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii)
372 
373            jsub=3*(iatoms(2,ii)-1)+jj
374            ksub=3*(iatoms(1,ii)-1)+kk
375            matrix(jsub,ksub)=matrix(jsub,ksub)-&
376 &           badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii)
377 
378 !          The diagonal blocks
379            jsub=3*(iatoms(1,ii)-1)+jj
380            ksub=3*(iatoms(1,ii)-1)+kk
381            matrix(jsub,ksub)=matrix(jsub,ksub)+&
382 &           badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii)
383 
384            jsub=3*(iatoms(2,ii)-1)+jj
385            ksub=3*(iatoms(2,ii)-1)+kk
386            matrix(jsub,ksub)=matrix(jsub,ksub)+&
387 &           badgerfactor*bonds%bond_vect(jj,ii)*bonds%bond_vect(kk,ii)
388 
389          end do !do kk=1,3
390        end do !do jj=1,3
391 
392      end if
393 
394    end do
395 
396    if (bonds%nbonds>0)then
397      ABI_DEALLOCATE(periods)
398      ABI_DEALLOCATE(iatoms)
399    end if
400 
401    call bonds_free(bonds)
402 
403    if (3*ab_mover%natom<100)then
404 !    Visualize the matrix
405      do jj=1,3*ab_mover%natom
406        write (std_out,fmt) matrix(jj,:)
407      end do
408    end if
409 
410    ABI_ALLOCATE(matrix_tmp,(3*ab_mover%natom,3*ab_mover%natom))
411    matrix_tmp(:,:)=matrix(:,:)
412    !write(*,*)"matrix_tmp",matrix_tmp
413 
414    ABI_ALLOCATE(work,(1))
415    lwork=-1
416    call DSYEV('V', 'U', 3*ab_mover%natom, matrix_tmp, 3*ab_mover%natom, w , work, lwork, info )
417    lwork=work(1)
418    write(std_out,*) '[DSYEV] Recommended lwork=',lwork
419    ABI_DEALLOCATE(work)
420    ABI_ALLOCATE(work,(lwork))
421    call DSYEV('V', 'U', 3*ab_mover%natom, matrix_tmp, 3*ab_mover%natom, w , work, lwork, info )
422    ABI_DEALLOCATE(work)
423    ABI_DEALLOCATE(matrix_tmp)
424    write(std_out,*) 'DSYEV info=',info
425    write(std_out,*) 'Eigenvalues:'
426    write(std_out,fmt) w(:)
427 
428    sigma=0
429    do jj=1,3*ab_mover%natom
430      sigma=max(w(jj),sigma)
431    end do
432 
433    matrix=lambda*matrix
434 
435    write(std_out,*) ch10
436    do ii=1,3*ab_mover%natom
437      matrix(ii,ii)=matrix(ii,ii)+(1-lambda)*sigma
438    end do
439 
440  end if ! if (Compute_Matrix)
441 
442 !##########################################################
443 !### 05. Use the precondition matrix to compute new residuals
444 
445  B=reshape(fcart,(/ 3*ab_mover%natom /))
446 
447  if (3*ab_mover%natom<100)then
448 !  Visualize the matrix
449    do jj=1,3*ab_mover%natom
450      write (std_out,fmt) matrix(jj,:)
451    end do
452  end if
453 
454 !call dsysv( uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info )
455 !MGNAG FIXME: This call causes a floating point exception if NAG+MKL
456  ABI_ALLOCATE(work,(1))
457  lwork=-1
458  call DSYSV( 'U', 3*ab_mover%natom, 1, matrix,&
459 & 3*ab_mover%natom, ipiv, B, 3*ab_mover%natom, work, lwork, info )
460 
461  lwork=work(1)
462  write(std_out,*) '[DSYSV] Recomended lwork=',lwork
463  ABI_DEALLOCATE(work)
464  ABI_ALLOCATE(work,(lwork))
465  call DSYSV( 'U', 3*ab_mover%natom, 1, matrix,&
466 & 3*ab_mover%natom, ipiv, B, 3*ab_mover%natom, work, lwork, info )
467  ABI_DEALLOCATE(work)
468 
469  write(std_out,*) 'DSYSV info=',info
470  write(std_out,*) 'Solution:'
471  write(std_out,fmt) B(:)
472 
473  forstr%fcart=reshape(B,(/ 3, ab_mover%natom /) )
474  call fcart2fred(forstr%fcart,forstr%fred,rprimd,ab_mover%natom)
475 
476 end subroutine prec_simple

ABINIT/pred_simple [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_simple

FUNCTION

 Ionmov predictors (4 & 6) Internal to SCFV

 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

PARENTS

      mover

CHILDREN

SOURCE

 81 subroutine pred_simple(ab_mover,hist,iexit)
 82 
 83 
 84 !This section has been created automatically by the script Abilint (TD).
 85 !Do not modify the following lines by hand.
 86 #undef ABI_FUNC
 87 #define ABI_FUNC 'pred_simple'
 88 !End of the abilint section
 89 
 90  implicit none
 91 
 92 !Arguments ------------------------------------
 93 !scalars
 94  type(abimover),intent(in) :: ab_mover
 95  type(abihist),intent(inout) :: hist
 96  integer,intent(in) :: iexit
 97 
 98 !Local variables-------------------------------
 99 !scalars
100  integer  :: ihist_next,kk,jj
101 
102 !***************************************************************************
103 !Beginning of executable session
104 !***************************************************************************
105 
106  if(iexit/=0)then
107    return
108  end if
109 
110 !All the operations are internal to scfcv.F90
111 
112 !XRED, FCART and VEL
113  ihist_next = abihist_findIndex(hist,+1)
114  do kk=1,ab_mover%natom
115    do jj=1,3
116      hist%xred(jj,kk, ihist_next)=hist%xred (jj,kk,hist%ihist)
117      hist%fcart(jj,kk,ihist_next)=hist%fcart(jj,kk,hist%ihist)
118      hist%vel(jj,kk,ihist_next)=hist%vel(jj,kk,hist%ihist)
119    end do
120  end do
121 
122 !ACELL
123  do jj=1,3
124    hist%acell(jj,ihist_next)=hist%acell(jj,hist%ihist)
125  end do
126 
127 !RPRIMD
128  do kk=1,3
129    do jj=1,3
130      hist%rprimd(jj,kk,ihist_next)=hist%rprimd(jj,kk,hist%ihist)
131    end do
132  end do
133 
134  hist%ihist=ihist_next
135 
136 end subroutine pred_simple