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 : !!***
|
