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