Line data Source code
1 : !!****m* ABINIT/m_pred_delocint
2 : !! NAME
3 : !! m_pred_delocint
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (MVer, DCA, XG, GMR, JCC, SE)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_pred_delocint
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 : use m_abimover
27 : use m_abihist
28 : use m_xfpack
29 : use m_linalg_interfaces
30 :
31 : use m_geometry, only : fcart2gred, xcart2xred, xred2xcart, metric, acrossb
32 : use m_bfgs, only : hessinit, hessupdt, brdene
33 : use m_results_gs, only : results_gs_type
34 :
35 : implicit none
36 :
37 : private
38 : !!***
39 :
40 : public :: pred_delocint
41 :
42 : contains
43 : !!***
44 :
45 : !!****f* ABINIT/pred_delocint
46 : !! NAME
47 : !! pred_delocint
48 : !!
49 : !! FUNCTION
50 : !! Ionmov predictors (10) BFGS with delocalized internal coordinates
51 : !!
52 : !! IONMOV 10:
53 : !! Given a starting point xred that is a vector of length 3*(natom-1)
54 : !! (reduced nuclei coordinates),
55 : !! and unit cell parameters (acell and rprimd) the
56 : !! Broyden-Fletcher-Goldfarb-Shanno minimization is performed on the
57 : !! total energy function, using its gradient (atomic forces and stresses)
58 : ! as calculated by the routine scfcv. Some atoms can be kept fixed,
59 : !! while the optimization of unit cell
60 : !! parameters is only performed if optcell/=0.
61 : !! The convergence requirement on
62 : !! the atomic forces, 'tolmxf', allows an early exit.
63 : !! Otherwise no more than 'ntime' steps are performed.
64 : !! Returned quantities are xred, and eventually acell and rprimd (new ones!).
65 : !! Could see Numerical Recipes (Fortran), 1986, page 307.
66 : !!
67 : !! Implements the delocalized internal coordinate scheme
68 : !! of Andzelm et al. in CPL .335. 321 (2001) \
69 : !! and Baker et al. JCP .105. 192 (1996)
70 : !!
71 : !! B matrix is derivative of delocalized internals wrt cartesian coordinates
72 : !! U matrix is eigenvectors of G = B*B^{T}
73 : !! S matrix is eigenvectors of F = B^{T}B
74 : !!
75 : !! INPUTS
76 : !! ab_mover <type(abimover)> : Datatype with all the information needed by the preditor
77 : !! itime : Index of the present iteration
78 : !! ntime : Maximal number of iterations
79 : !! ionmov : (10 or 11) Specific kind of BFGS
80 : !! zDEBUG : if true print some debugging information
81 : !!
82 : !! OUTPUT
83 : !!
84 : !! SIDE EFFECTS
85 : !! hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
86 : !!
87 : !! SOURCE
88 :
89 2 : subroutine pred_delocint(ab_mover,ab_xfh,deloc,forstr,hist,ionmov,itime,zDEBUG,iexit)
90 :
91 : !Arguments ------------------------------------
92 : !scalars
93 : type(abimover),intent(in) :: ab_mover
94 : type(ab_xfh_type),intent(inout) :: ab_xfh
95 : type(abihist),intent(inout) :: hist
96 : type(abiforstr),intent(in) :: forstr
97 : type(delocint),intent(inout) :: deloc
98 : integer,intent(in) :: itime
99 : integer,intent(in) :: ionmov
100 : integer,intent(in) :: iexit
101 : logical,intent(in) :: zDEBUG
102 :
103 : !Local variables-------------------------------
104 : !scalars
105 : integer :: ndim,cycl_main
106 : integer :: ihist_prev,ii,jj,kk
107 : real(dp),save :: ucvol0
108 : real(dp) :: ucvol
109 : real(dp) :: etotal,etotal_prev
110 : logical :: DEBUG=.TRUE.
111 : !integer,save :: icenter,irshift ! DELOCINT indexes
112 : integer,save :: ndeloc ! DELOCINT number of
113 : character(len=500) :: message
114 :
115 : !arrays
116 : real(dp),allocatable,save :: hessin(:,:),vin(:),vin_prev(:)
117 : real(dp),allocatable,save :: vout(:),vout_prev(:)
118 : real(dp),save :: acell0(3) ! Initial acell
119 : real(dp),save :: rprimd0(3,3) ! Initial rprimd
120 2 : real(dp),allocatable :: prim_int(:)
121 : real(dp),allocatable,save :: u_matrix(:,:) ! DELOCINT this may need to be added to type inside ab_mover
122 : real(dp) :: acell(3)
123 : real(dp) :: rprimd(3,3)
124 : real(dp) :: gprimd(3,3)
125 : real(dp) :: gmet(3,3)
126 : real(dp) :: rmet(3,3)
127 4 : real(dp) :: residual(3,ab_mover%natom)
128 : !real(dp) :: residual_corrected(3,ab_mover%natom)
129 4 : real(dp) :: xred(3,ab_mover%natom),xcart(3,ab_mover%natom)
130 : real(dp) :: strten(6)
131 4 : real(dp) :: deloc_gred(3*(ab_mover%natom-1))
132 4 : real(dp) :: deloc_int(3*(ab_mover%natom-1))
133 2 : real(dp) :: bt_inv_matrix(3*(ab_mover%natom-1),3*ab_mover%natom)
134 :
135 : !***************************************************************************
136 : !Beginning of executable session
137 : !***************************************************************************
138 :
139 2 : if(iexit/=0)then
140 1 : ABI_SFREE(vin)
141 1 : ABI_SFREE(vout)
142 1 : ABI_SFREE(vin_prev)
143 1 : ABI_SFREE(vout_prev)
144 1 : ABI_SFREE(hessin)
145 1 : ABI_SFREE(u_matrix)
146 : return
147 : end if
148 :
149 : !write(std_out,*) 'delocint 01'
150 : !##########################################################
151 : !### 01. Debugging and Verbose
152 :
153 1 : if(DEBUG)then
154 43 : write(std_out,'(a,3a,38a,39a)') ch10,('-',kk=1,3),'Debugging and Verbose for pred_deloint',('-',kk=1,39)
155 1 : write(std_out,*) 'ionmov: ',ionmov
156 1 : write(std_out,*) 'itime: ',itime
157 : end if
158 :
159 : !write(std_out,*) 'delocint 02'
160 : !##########################################################
161 : !### 02. Compute the dimension of vectors (ndim)
162 :
163 : !With internal we have 1 coordinate less
164 1 : ndeloc = 3*(ab_mover%natom-1)
165 1 : ndim=ndeloc
166 4 : deloc_int(:)=zero
167 4 : deloc_gred(:)=zero
168 1 : if(ab_mover%optcell==1) ndim=ndim+1
169 1 : if(ab_mover%optcell==2 .or.&
170 0 : & ab_mover%optcell==3) ndim=ndim+6
171 1 : if(ab_mover%optcell>=4) ndim=ndim+3
172 :
173 1 : if(DEBUG) write(std_out,*) 'Dimension of vin, vout and hessian (ndim): ',ndim
174 :
175 : !write(std_out,*) 'delocint 03'
176 : !##########################################################
177 : !### 03. Allocate the vectors vin, vout and hessian matrix
178 :
179 : !Notice thqt vin, vout, etc could be allocated
180 : !From a previous dataset with a different ndim
181 1 : if(itime==1)then
182 1 : ABI_SFREE(vin)
183 1 : ABI_SFREE(vout)
184 1 : ABI_SFREE(vin_prev)
185 1 : ABI_SFREE(vout_prev)
186 1 : ABI_SFREE(hessin)
187 :
188 3 : ABI_MALLOC(vin,(ndim))
189 2 : ABI_MALLOC(vout,(ndim))
190 2 : ABI_MALLOC(vin_prev,(ndim))
191 2 : ABI_MALLOC(vout_prev,(ndim))
192 4 : ABI_MALLOC(hessin,(ndim,ndim))
193 : end if
194 :
195 :
196 : !write(std_out,*) 'delocint 04'
197 : !##########################################################
198 : !### 04. Obtain the present values from the history
199 :
200 1 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
201 1 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
202 :
203 7 : strten(:)=hist%strten(:,hist%ihist)
204 1 : etotal =hist%etot(hist%ihist)
205 :
206 : !Fill the residual with forces (No preconditioning)
207 : !Or the preconditioned forces
208 1 : if (ab_mover%goprecon==0)then
209 1 : call fcart2gred(hist%fcart(:,:,hist%ihist),residual,rprimd,ab_mover%natom)
210 : else
211 0 : residual(:,:)= forstr%gred(:,:)
212 : end if
213 :
214 1 : if(zDEBUG)then
215 0 : write (std_out,*) 'residual:'
216 0 : do kk=1,ab_mover%natom
217 0 : write (std_out,*) residual(:,kk)
218 : end do
219 0 : write (std_out,*) 'strten:'
220 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
221 0 : write (std_out,*) 'etotal:'
222 0 : write (std_out,*) etotal
223 : end if
224 :
225 1 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
226 :
227 : !Save initial values
228 1 : if (itime==1)then
229 1 : acell0(:)=acell(:)
230 1 : rprimd0(:,:)=rprimd(:,:)
231 1 : ucvol0=ucvol
232 : end if
233 :
234 : !DEBUG (UCVOL)
235 1 : if(DEBUG)then
236 1 : write(std_out,*) 'Volume of cell (ucvol):',ucvol
237 : end if
238 :
239 : !Get rid of mean force on whole unit cell, but only if no
240 : !generalized constraints are in effect
241 : ! residual_corrected(:,:)=residual(:,:)
242 : ! if(ab_mover%nconeq==0)then
243 : ! do ii=1,3
244 : ! if (ii/=3.or.ab_mover%jellslab==0) then
245 : ! favg=sum(residual_corrected(ii,:))/dble(ab_mover%natom)
246 : ! residual_corrected(ii,:)=residual_corrected(ii,:)-favg
247 : ! end if
248 : ! end do
249 : ! end if
250 :
251 : !write(std_out,*) 'delocint 05'
252 : !##########################################################
253 : !### 05. Compute internals for first time
254 :
255 1 : if (itime==1)then
256 : call make_prim_internals(deloc,ab_mover%natom,&
257 1 : & ab_mover%ntypat,rprimd,ab_mover%typat,xcart,ab_mover%znucl)
258 :
259 3 : ABI_MALLOC(prim_int,(deloc%ninternal))
260 :
261 1 : if(DEBUG)then
262 1 : write (message,'(a,i6)') 'Number of primitive internal coordinates (ninternal): ',deloc%ninternal
263 1 : call wrtout(std_out, message,'COLL')
264 : end if
265 :
266 1 : if (allocated(u_matrix)) then
267 0 : ABI_FREE(u_matrix)
268 : end if
269 4 : ABI_MALLOC(u_matrix,(deloc%ninternal,ndeloc))
270 :
271 1 : call calc_prim_int(deloc,ab_mover%natom,rprimd,xcart,prim_int)
272 :
273 1 : if(DEBUG)then
274 1 : write (message,'(a)') 'Primitive internal coordinate values:'
275 1 : call wrtout(std_out, message,'COLL')
276 1 : write (message,'(a)') ' Bonds:'
277 1 : call wrtout(std_out, message,'COLL')
278 5 : do ii = 1, deloc%nbond
279 4 : write (message,'(i6,E20.10)') ii, prim_int(ii)
280 5 : call wrtout(std_out, message,'COLL')
281 : end do
282 :
283 1 : write (message,'(a)') ' Angles:'
284 1 : call wrtout(std_out, message,'COLL')
285 13 : do ii = deloc%nbond+1, deloc%nbond+deloc%nang
286 12 : write (message,'(i6,2(E20.10,2x))') ii, prim_int(ii), prim_int(ii)/pi*180.0_dp
287 13 : call wrtout(std_out, message,'COLL')
288 : end do
289 :
290 1 : write (message,'(a)') ' Dihedrals:'
291 1 : call wrtout(std_out, message,'COLL')
292 37 : do ii = deloc%nbond+deloc%nang+1, deloc%nbond+deloc%nang+deloc%ndihed
293 36 : write (message,'(i6,2(E20.10,2x))') ii, prim_int(ii), prim_int(ii)/pi*180.0_dp
294 37 : call wrtout(std_out, message,'COLL')
295 : end do
296 :
297 1 : write (message,'(a)') ' Cartesian auxiliary coordinates for constraints:'
298 1 : call wrtout(std_out, message,'COLL')
299 1 : do ii = deloc%nbond+deloc%nang+deloc%ndihed+1, deloc%ninternal
300 0 : write (message,'(i6,E20.10)') ii, prim_int(ii)
301 1 : call wrtout(std_out, message,'COLL')
302 : end do
303 : end if
304 :
305 1 : ABI_FREE(prim_int)
306 :
307 : ! equal weight on all internal coordinates as a starting point.
308 160 : u_matrix(:,:) = one / dble (ndeloc)
309 :
310 : ! Zero the arrays before first use
311 4 : deloc_gred(:) = zero
312 :
313 : end if
314 :
315 3 : ABI_MALLOC(prim_int,(deloc%ninternal))
316 :
317 : !write(std_out,*) 'delocint 06'
318 : !##########################################################
319 : !### 06. Compute delocalized coordinates and forces
320 :
321 : !xcart ---> deloc_int
322 :
323 : !Convert positions to delocalized coordinates for next step
324 : call xcart2deloc(deloc,ab_mover%natom,rprimd,xcart,&
325 1 : & bt_inv_matrix,u_matrix,deloc_int,prim_int)
326 :
327 : !gred ---> deloc_gred
328 :
329 : !Convert gradients to delocalized coordinates for next step
330 1 : call gred2gdeloc(bt_inv_matrix,deloc_gred,residual,ab_mover%natom,gprimd)
331 :
332 : !write(std_out,*) 'delocint 07'
333 : !##########################################################
334 : !### 07. Fill the vectors vin and vout
335 :
336 : !DEBUG deloc_int and deloc_gred before pack
337 1 : if(DEBUG)then
338 1 : write (std_out,*) 'Delocalized internals and forces (ndeloc):',ndeloc
339 1 : write(std_out,*) 'deloc_int'
340 2 : do ii=1,ndeloc,3
341 2 : if (ii+2<=ndeloc)then
342 1 : write(std_out,*) ii,deloc_int(ii:ii+2)
343 : else
344 0 : write(std_out,*) ii,deloc_int(ii:ndeloc)
345 : end if
346 : end do
347 1 : write(std_out,*) 'deloc_gred'
348 2 : do ii=1,ndeloc,3
349 2 : if (ii+2<=ndeloc)then
350 1 : write(std_out,*) ii,deloc_gred(ii:ii+2)
351 : else
352 0 : write(std_out,*) ii,deloc_gred(ii:ndeloc)
353 : end if
354 : end do
355 : end if
356 :
357 : !DELOCINT
358 : !Instead of gred_corrected we use deloc_gred
359 : !Instead of xred e use deloc_int
360 : !
361 : !Initialize input vectors : first vin, then vout
362 : !The values of vin from the previous iteration
363 : !should be the same
364 : call xfpack_x2vin(acell, ab_mover%natom-1, ndim,&
365 : & ab_mover%nsym, ab_mover%optcell, rprimd, rprimd0,&
366 1 : & ab_mover%symrel, ucvol, ucvol0, vin, deloc_int)
367 : !end if
368 :
369 : call xfpack_f2vout(deloc_gred, ab_mover%natom-1, ndim,&
370 : & ab_mover%optcell, ab_mover%strtarget, strten, ucvol, &
371 1 : & vout)
372 :
373 : !write(std_out,*) 'delocint 08'
374 : !##########################################################
375 : !### 08. Initialize or update the hessian matrix
376 :
377 : !Initialise the Hessian matrix using gmet
378 1 : if (itime==1)then
379 :
380 : ! Initialise the Hessian matrix with ab_mover%userrc.
381 : ! this has become unusable because it imposes ndim >= 3 natom
382 : ! ident = 3x3 identity matrix
383 : ! call hessinit(ab_mover, hessin, gmet, ndim, ucvol)
384 13 : hessin = zero
385 4 : do ii=1, ndim
386 4 : hessin (ii,ii) = one
387 : end do
388 :
389 : ! ! Initialize inverse hessian with identity matrix
390 : ! ! in cartesian coordinates, which makes use of metric tensor gmet
391 : ! ! in reduced coordinates.
392 : ! hessin(:,:)=zero
393 : ! do ii=1,ab_mover%natom
394 : ! do kk=1,3
395 : ! do jj=1,3
396 : ! ! Warning : implemented in reduced coordinates
397 : ! if (ab_mover%iatfix(kk,ii)==0 .and.&
398 : ! & ab_mover%iatfix(jj,ii)==0 )then
399 : ! hessin(kk+3*(ii-1),jj+3*(ii-1))=gmet(kk,jj)
400 : ! end if
401 : ! end do
402 : ! end do
403 : ! end do
404 : ! if(ab_mover%optcell/=0)then
405 : ! ! These values might lead to too large changes in some cases
406 : ! diag=ab_mover%strprecon*30.0_dp/ucvol
407 : ! if(ab_mover%optcell==1) diag=diag/three
408 : ! do ii=3*ab_mover%natom+1,ndim
409 : ! hessin(ii,ii)=diag
410 : ! end do
411 : ! end if
412 :
413 1 : if (ab_mover%restartxf/=0) then
414 :
415 : call xfh_recover_deloc(ab_xfh,ab_mover,acell,cycl_main,&
416 : & residual,hessin,ndim,rprimd,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,&
417 : & vout,vout_prev,xred,deloc,deloc_int,deloc_gred,bt_inv_matrix,gprimd,prim_int,&
418 0 : & u_matrix)
419 :
420 : end if
421 :
422 : end if
423 :
424 1 : ABI_FREE(prim_int)
425 :
426 1 : if(itime>1)then
427 : ! Update the hessian matrix, by taking into account the
428 : ! current pair (x,f) and the previous one.
429 : call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,vin,&
430 0 : & vin_prev,vout,vout_prev)
431 :
432 : end if
433 :
434 : !DEBUG (vin,vout and hessin before prediction)
435 1 : if(DEBUG)then
436 1 : write(std_out,*) 'Vectors vin and vout and inverse of Hessian (hessin) [before]'
437 1 : write(std_out,*) 'vin:'
438 2 : do ii=1,ndim,3
439 2 : if (ii+2<=ndim)then
440 1 : write(std_out,*) ii,vin(ii:ii+2)
441 : else
442 0 : write(std_out,*) ii,vin(ii:ndim)
443 : end if
444 : end do
445 1 : write(std_out,*) 'vout:'
446 2 : do ii=1,ndim,3
447 2 : if (ii+2<=ndim)then
448 1 : write(std_out,*) ii,vout(ii:ii+2)
449 : else
450 0 : write(std_out,*) ii,vout(ii:ndim)
451 : end if
452 : end do
453 1 : write(std_out,*) 'Inverse Hessian (hessin): ',ndim,'x',ndim
454 4 : do kk=1,ndim
455 7 : do jj=1,ndim,3
456 6 : if (jj+2<=ndim)then
457 3 : write(std_out,*) jj,hessin(jj:jj+2,kk)
458 : else
459 0 : write(std_out,*) jj,hessin(jj:ndim,kk)
460 : end if
461 : end do
462 : end do
463 : end if
464 :
465 : !write(std_out,*) 'delocint 09'
466 : !##########################################################
467 : !### 09. Compute the next values
468 :
469 1 : if(ionmov==10 .or. itime==1)then
470 :
471 : ! Previous cartesian coordinates
472 4 : vin_prev(:)=vin(:)
473 :
474 : ! New atomic cartesian coordinates are obtained from vin, hessin
475 : ! and vout
476 40 : vin = vin - matmul(hessin, vout)
477 : ! Previous atomic forces
478 4 : vout_prev(:)=vout(:)
479 :
480 : else
481 0 : if(ionmov==11)then
482 0 : ihist_prev = abihist_findIndex(hist,-1)
483 0 : etotal_prev=hist%etot(ihist_prev)
484 : ! Here the BFGS algorithm, modified to take into account the
485 : ! energy
486 : call brdene(etotal,etotal_prev,hessin,&
487 0 : & ndim,vin,vin_prev,vout,vout_prev)
488 :
489 : end if
490 :
491 : ! DEBUG (vin,vout and hessin after prediction)
492 0 : if(DEBUG)then
493 0 : write(std_out,*) 'Vectors vin and vout [after prediction]'
494 0 : write(std_out,*) 'vin:'
495 0 : do ii=1,ndim,3
496 0 : if (ii+2<=ndim)then
497 0 : write(std_out,*) ii,vin(ii:ii+2)
498 : else
499 0 : write(std_out,*) ii,vin(ii:ndim)
500 : end if
501 : end do
502 0 : write(std_out,*) 'vout:'
503 0 : do ii=1,ndim,3
504 0 : if (ii+2<=ndim)then
505 0 : write(std_out,*) ii,vout(ii:ii+2)
506 : else
507 0 : write(std_out,*) ii,vout(ii:ndim)
508 : end if
509 : end do
510 : end if
511 :
512 : ! Implement fixing of atoms : put back old values for fixed
513 : ! components
514 0 : do kk=1,ab_mover%natom
515 0 : do jj=1,3
516 : ! Warning : implemented in reduced coordinates
517 0 : if ( ab_mover%iatfix(jj,kk)==1) then
518 0 : vin(jj+(kk-1)*3)=vin_prev(jj+(kk-1)*3)
519 : end if
520 : end do
521 : end do
522 : end if
523 :
524 : !write(std_out,*) 'delocint 10'
525 : !##########################################################
526 : !### 10. Convert from delocalized to xcart and xred
527 :
528 : !Transfer vin to deloc_int, acell and rprimd
529 : call xfpack_vin2x(acell, acell0, ab_mover%natom-1, ndim,&
530 : & ab_mover%nsym, ab_mover%optcell, rprimd, rprimd0,&
531 : & ab_mover%symrel, ucvol, ucvol0,&
532 1 : & vin, deloc_int)
533 :
534 1 : if(DEBUG)then
535 1 : write (std_out,*) 'Delocalized internals (deloc_int) [after prediction]:'
536 1 : write(std_out,*) 'deloc_int:'
537 2 : do ii=1,ndeloc,3
538 2 : if (ii+2<=ndeloc)then
539 1 : write(std_out,*) ii,deloc_int(ii:ii+2)
540 : else
541 0 : write(std_out,*) ii,deloc_int(ii:ndeloc)
542 : end if
543 : end do
544 1 : write(std_out,*) 'BT Inverse Matrix:'
545 7 : do ii=1,3*ab_mover%natom
546 7 : write(std_out,*) bt_inv_matrix(:,ii)
547 : end do
548 1 : write (std_out,*) 'xcart (before deloc2xcart):'
549 3 : do ii=1,ab_mover%natom
550 3 : write (std_out,*) xcart(:,ii)
551 : end do
552 : end if
553 :
554 : !this routine contains an iterative scheme to find xcart
555 : !from the non-linear relations between deloc and xcart
556 : !SIGNIFICANTLY DIFFERENT FROM xcart2deloc
557 : call deloc2xcart(deloc,ab_mover%natom,rprimd,xcart,&
558 1 : & deloc_int,bt_inv_matrix,u_matrix)
559 :
560 : !Convert new xcart (cartesian) to xred (reduced coordinates)
561 1 : call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
562 :
563 :
564 : !write(std_out,*) 'delocint 11'
565 : !##########################################################
566 : !### 11. Update the history with the prediction
567 :
568 : !Increase indexes
569 1 : hist%ihist = abihist_findIndex(hist,+1)
570 :
571 1 : if(ab_mover%optcell/=0)then
572 0 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
573 : end if
574 :
575 : !Fill the history with the variables
576 : !xred, acell, rprimd, vel
577 1 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
578 1 : ihist_prev = abihist_findIndex(hist,-1)
579 9 : hist%vel(:,:,hist%ihist)=hist%vel(:,:,ihist_prev)
580 :
581 1 : if(zDEBUG)then
582 0 : write (std_out,*) 'residual:'
583 0 : do kk=1,ab_mover%natom
584 0 : write (std_out,*) residual(:,kk)
585 : end do
586 0 : write (std_out,*) 'strten:'
587 0 : write (std_out,*) strten(1:3),ch10,strten(4:6)
588 0 : write (std_out,*) 'etotal:'
589 0 : write (std_out,*) etotal
590 : end if
591 :
592 : end subroutine pred_delocint
593 : !!***
594 :
595 : !!****f* ABINIT/deloc2xcart
596 : !! NAME
597 : !! deloc2xcart
598 : !!
599 : !! FUNCTION
600 : !! Determine the cartesian coordinates which correspond to the
601 : !! given values of the delocalized coordinates. The relationship
602 : !! is non-linear, so use an iterative scheme, as in Baker
603 : !! JCP .105. 192 (1996).
604 : !! Older reference: Pulay and co. JACS 101 2550 (1979)
605 : !!
606 : !! INPUTS
607 : !! deloc <type(delocint)>=Important variables for
608 : !! | pred_delocint
609 : !! |
610 : !! | nang = Number of angles
611 : !! | nbond = Number of bonds
612 : !! | ncart = Number of cartesian directions
613 : !! | (used for constraints)
614 : !! | ndihed = Number of dihedrals
615 : !! | nrshift = Dimension of rshift
616 : !! | ninternal= Number of internal coordinates
617 : !! | ninternal=nbond+nang+ndihed+ncart
618 : !! |
619 : !! | angs(2,3,nang) = Indexes to characterize angles
620 : !! | bonds(2,2,nbond)= For a bond between iatom and jatom
621 : !! | bonds(1,1,nbond) = iatom
622 : !! | bonds(2,1,nbond) = icenter
623 : !! | bonds(1,2,nbond) = jatom
624 : !! | bonds(2,2,nbond) = irshift
625 : !! | carts(2,ncart) = Index of total primitive internal,
626 : !! | and atom (carts(2,:))
627 : !! | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals
628 : !! |
629 : !! | rshift(3,nrshift)= Shift in xred that must be done to find
630 : !! | all neighbors of a given atom within a
631 : !! | given number of neighboring shells
632 : !! natom = Number of atoms (dtset%natom)
633 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
634 : !!
635 : !! OUTPUT
636 : !! bt_inv_matrix(3*(natom-1),3*natom)=inverse of transpose of B matrix
637 : !!
638 : !! SIDE EFFECTS
639 : !! u_matrix(ninternal,3*(natom-1))=eigenvectors of G = BB^T matrix
640 : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
641 : !!
642 : !! NOTES
643 : !!
644 : !! SOURCE
645 :
646 1 : subroutine deloc2xcart(deloc,natom,rprimd,xcart,deloc_int,btinv,u_matrix)
647 :
648 : !Arguments ------------------------------------
649 : !scalars
650 : integer,intent(in) :: natom
651 : type(delocint),intent(in) :: deloc
652 : !arrays
653 : real(dp),intent(in) :: deloc_int(3*(natom-1)),rprimd(3,3)
654 : real(dp),intent(inout) :: u_matrix(deloc%ninternal,3*(natom-1))
655 : real(dp),intent(inout) :: xcart(3,natom)
656 : real(dp),intent(out) :: btinv(3*(natom-1),3*natom)
657 :
658 : !Local variables-------------------------------
659 : !scalars
660 : integer :: iiter,iprim,niter
661 : integer :: ii
662 : real(dp) :: minmix, maxmix
663 : real(dp) :: mix,tot_diff, toldeloc
664 : real(dp) :: lntoldeloc
665 : logical :: DEBUG=.FALSE.
666 : !arrays
667 2 : real(dp) :: btinv_tmp(3*(natom-1),3*natom)
668 2 : real(dp) :: cgrad(3*natom),cgrad_old(3*natom)
669 2 : real(dp) :: deloc_int_now(3*(natom-1)),prim_int(deloc%ninternal)
670 2 : real(dp) :: tmpxcart(3*natom)
671 2 : real(dp) :: xdeloc_diff(3*(natom-1))
672 :
673 : character(len=500) :: message
674 :
675 : ! ******************************************************************
676 :
677 1 : if (DEBUG) then
678 0 : write(ab_out,*) 'ENTERING DELOC2XCART'
679 :
680 0 : write (message,*) 'BONDS=',deloc%nbond
681 0 : call wrtout(ab_out,message,'COLL')
682 0 : do ii = 1, deloc%nbond
683 0 : write (message,*) ii, deloc%bonds(:,:,ii)
684 0 : call wrtout(ab_out,message,'COLL')
685 : end do
686 :
687 0 : write (message,*) 'ANGS=',deloc%nang
688 0 : call wrtout(ab_out,message,'COLL')
689 0 : do ii = 1, deloc%nang
690 0 : write (message,*) ii, deloc%angs(:,:,ii)
691 0 : call wrtout(ab_out,message,'COLL')
692 : end do
693 :
694 0 : write (message,*) 'DIHEDRALS=',deloc%ndihed
695 0 : call wrtout(ab_out,message,'COLL')
696 0 : do ii = 1, deloc%ndihed
697 0 : write (message,*) ii, deloc%dihedrals(:,:,ii)
698 0 : call wrtout(ab_out,message,'COLL')
699 : end do
700 :
701 0 : write (message,*) 'CARTS=',deloc%ncart
702 0 : call wrtout(ab_out,message,'COLL')
703 0 : do ii = 1, deloc%ncart
704 0 : write (message,*) ii, deloc%carts(:,ii)
705 0 : call wrtout(ab_out,message,'COLL')
706 : end do
707 :
708 0 : write (ab_out,*) 'xcart (input)'
709 0 : do ii=1,natom
710 0 : write (ab_out,*) xcart(:,ii)
711 : end do
712 :
713 : end if
714 :
715 1 : niter = 200
716 2 : tmpxcart = reshape(xcart,(/3*natom/))
717 :
718 7 : cgrad_old(:) = zero
719 7 : cgrad(:) = zero
720 : maxmix = 0.9_dp
721 : minmix = 0.2_dp
722 : toldeloc = tol10
723 : lntoldeloc = log(toldeloc)
724 :
725 62 : do iiter=1,niter
726 62 : if (iiter==1) then
727 : mix= minmix
728 : else
729 61 : mix = minmix + (maxmix-minmix)*(log(tot_diff)-lntoldeloc) / lntoldeloc
730 : end if
731 61 : if (mix < minmix) mix = minmix
732 : if (mix > maxmix) mix = maxmix
733 :
734 434 : tmpxcart(:) = tmpxcart(:) + mix*cgrad(:)
735 186 : xcart = reshape(tmpxcart,(/3,natom/))
736 : call xcart2deloc(deloc,natom,rprimd,xcart,&
737 62 : & btinv_tmp,u_matrix,deloc_int_now,prim_int)
738 : ! update the BT^{-1} matrix?
739 1550 : btinv(:,:) = btinv_tmp(:,:)
740 :
741 248 : xdeloc_diff(:) = deloc_int(:) - deloc_int_now(:)
742 :
743 248 : tot_diff = sum(abs(xdeloc_diff))
744 62 : if (tot_diff < toldeloc) exit
745 :
746 427 : cgrad_old(:) = cgrad(:)
747 :
748 : ! gradient vector = btinv^{T} * xdeloc_diff
749 : call dgemv('T',3*(natom-1),3*natom,one,&
750 62 : & btinv,3*(natom-1),xdeloc_diff,1,zero,cgrad,1)
751 : end do
752 : !end iiter do
753 :
754 : call xcart2deloc(deloc,natom,rprimd,xcart,&
755 1 : & btinv,u_matrix,deloc_int_now,prim_int)
756 1 : write (message,'(3a)') 'delocalized internals, after convergence of xcart = ', ch10
757 1 : call wrtout(std_out,message,'COLL')
758 4 : do ii = 1, 3*(natom-1)
759 3 : write (message,'(I6,E20.10,2x)') ii, deloc_int_now(ii)
760 4 : call wrtout(std_out,message,'COLL')
761 : end do
762 :
763 4 : xdeloc_diff(:) = deloc_int(:) - deloc_int_now(:)
764 :
765 1 : write (message,'(a)') 'Primitive internal coordinate values:'
766 1 : call wrtout(std_out,message,'COLL')
767 5 : do iprim = 1, deloc%nbond
768 4 : write (message,'(i6,E20.10)') iprim, prim_int(iprim)
769 5 : call wrtout(std_out,message,'COLL')
770 : end do
771 49 : do iprim = deloc%nbond+1, deloc%nbond+deloc%nang+deloc%ndihed
772 48 : write (message,'(i6,2E20.10)') iprim, prim_int(iprim), prim_int(iprim)/pi*180.0_dp
773 49 : call wrtout(std_out,message,'COLL')
774 : end do
775 1 : do iprim = deloc%nbond+deloc%nang+deloc%ndihed+1, deloc%ninternal
776 0 : write (message,'(i6,E20.10)') iprim, prim_int(iprim)
777 1 : call wrtout(std_out,message,'COLL')
778 : end do
779 :
780 1 : if (iiter == niter+1) then
781 0 : write (message,'(a,i6,a,E20.10)') 'deloc2xcart : Error, xcart not converged in ', niter, 'iterations ', tot_diff
782 0 : ABI_ERROR(message)
783 : end if
784 :
785 1 : if(DEBUG)then
786 0 : write (ab_out,*) 'xcart (output)'
787 0 : do ii=1,natom
788 0 : write (ab_out,*) xcart(:,ii)
789 : end do
790 0 : write(ab_out,*) 'EXITING DELOC2XCART'
791 : end if
792 :
793 1 : end subroutine deloc2xcart
794 : !!***
795 :
796 : !!****f* ABINIT/gred2gdeloc
797 : !! NAME
798 : !! gred2gdeloc
799 : !!
800 : !! FUNCTION
801 : !! calculate delocalized forces from reduced coordinate ones
802 : !!
803 : !! INPUTS
804 : !! btinv(3*(natom-1),3*natom)= inverse transpose of B matrix (see delocint)
805 : !! natom = number of atoms
806 : !! gprimd(3,3)=dimensional translations in reciprocal space (bohr-1)
807 : !!
808 : !! OUTPUT
809 : !! deloc_gred(3*(natom-1))=delocalized gradients from reduced coordinate ones
810 : !! gred(3,natom)=delocalized gradients in reduced coordinates
811 : !!
812 : !! SOURCE
813 :
814 1 : subroutine gred2gdeloc(btinv,deloc_gred,gred,natom,gprimd)
815 :
816 : !Arguments ------------------------------------
817 : !scalars
818 : integer, intent(in) :: natom
819 : !arrays
820 : real(dp),intent(in) :: btinv(3*(natom-1),3*natom),gprimd(3,3),gred(3,natom)
821 : real(dp),intent(out) :: deloc_gred(3*(natom-1))
822 :
823 : !Local variables-------------------------------
824 : integer :: ii
825 : !arrays
826 2 : real(dp) :: fcart(3,natom)
827 : character(len=500) :: message
828 :
829 : ! ******************************************************************
830 :
831 : !make cartesian forces
832 :
833 : call dgemm('N','N',3,natom,3,one,&
834 1 : & gprimd,3,gred,3,zero,fcart,3)
835 :
836 : !turn cartesian to delocalized forces
837 : call dgemv('N',3*(natom-1),3*natom,one,&
838 1 : & btinv,3*(natom-1),fcart,1,zero,deloc_gred,1)
839 :
840 1 : write (message,'(a)') 'gred2gdeloc : deloc_gred = '
841 1 : call wrtout(std_out,message,'COLL')
842 :
843 4 : do ii = 1, 3*(natom-1)
844 3 : write (message,'(I6,E16.6)') ii, deloc_gred(ii)
845 4 : call wrtout(std_out,message,'COLL')
846 : end do
847 :
848 1 : end subroutine gred2gdeloc
849 : !!***
850 :
851 : !!****f* ABINIT/calc_b_matrix
852 : !! NAME
853 : !! calc_b_matrix
854 : !!
855 : !! FUNCTION
856 : !! calculate values of derivatives of internal coordinates as a function of
857 : !! cartesian ones = B matrix
858 : !!
859 : !! INPUTS
860 : !! angs= number of angles
861 : !! bonds(2,2,nbond)=for a bond between iatom and jatom
862 : !! bonds(1,1,nbond) = iatom
863 : !! bonds(2,1,nbond) = icenter
864 : !! bonds(1,2,nbond) = jatom
865 : !! bonds(2,2,nbond) = irshift
866 : !! carts(2,ncart)= index of total primitive internal, and atom (carts(2,:))
867 : !! dihedrals(2,4,ndihed)=indexes to characterize dihedrals
868 : !! nang(2,3,nang)=indexes to characterize angles
869 : !! nbond=number of bonds
870 : !! ncart=number of auxiliary cartesian atom coordinates (used for constraints)
871 : !! ndihed= number of dihedrals
872 : !! ninternal=nbond+nang+ndihed+ncart: number of internal coordinates
873 : !! nrshift= dimension of rshift
874 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
875 : !! rshift(3,nrshift)=shift in xred that must be done to find all neighbors of
876 : !! a given atom within a given number of neighboring shells
877 : !! xcart(3,natom)=cartesian coordinates of atoms (bohr)
878 : !!
879 : !! OUTPUT
880 : !! b_matrix(ninternal,3*natom)=matrix of derivatives of internal coordinates
881 : !! wrt cartesians
882 : !!
883 : !! SOURCE
884 :
885 64 : subroutine calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix)
886 :
887 : !Arguments ------------------------------------
888 : !scalars
889 : integer,intent(in) :: natom
890 : type(delocint),intent(in) :: deloc
891 :
892 : !arrays
893 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
894 : real(dp),intent(out) :: b_matrix(deloc%ninternal,3*natom)
895 :
896 : !Local variables-------------------------------
897 : !scalars
898 : integer :: i1,i2,i3,i4,iang,ibond,icart,idihed,iprim,s1,s2,s3,s4
899 : !arrays
900 : real(dp) :: bb(3),r1(3),r2(3),r3(3),r4(3)
901 :
902 : ! *************************************************************************
903 :
904 64 : iprim=0
905 20416 : b_matrix(:,:) = zero
906 :
907 320 : do ibond=1,deloc%nbond
908 256 : i1 = deloc%bonds(1,1,ibond)
909 256 : s1 = deloc%bonds(2,1,ibond)
910 : r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
911 : & +deloc%rshift(2,s1)*rprimd(:,2)&
912 1024 : & +deloc%rshift(3,s1)*rprimd(:,3)
913 256 : i2 = deloc%bonds(1,2,ibond)
914 256 : s2 = deloc%bonds(2,2,ibond)
915 : r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
916 : & +deloc%rshift(2,s2)*rprimd(:,2)&
917 1024 : & +deloc%rshift(3,s2)*rprimd(:,3)
918 256 : iprim=iprim+1
919 256 : call dbond_length_d1(r1,r2,bb)
920 1024 : b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
921 256 : call dbond_length_d1(r2,r1,bb)
922 1088 : b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
923 : end do
924 :
925 : !second: angle values (ang)
926 832 : do iang=1,deloc%nang
927 768 : i1 = deloc%angs(1,1,iang)
928 768 : s1 = deloc%angs(2,1,iang)
929 : r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
930 : & +deloc%rshift(2,s1)*rprimd(:,2)&
931 3072 : & +deloc%rshift(3,s1)*rprimd(:,3)
932 768 : i2 = deloc%angs(1,2,iang)
933 768 : s2 = deloc%angs(2,2,iang)
934 : r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
935 : & +deloc%rshift(2,s2)*rprimd(:,2)&
936 3072 : & +deloc%rshift(3,s2)*rprimd(:,3)
937 768 : i3 = deloc%angs(1,3,iang)
938 768 : s3 = deloc%angs(2,3,iang)
939 : r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
940 : & +deloc%rshift(2,s3)*rprimd(:,2)&
941 3072 : & +deloc%rshift(3,s3)*rprimd(:,3)
942 768 : iprim=iprim+1
943 768 : call dang_d1(r1,r2,r3,bb)
944 3072 : b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
945 768 : call dang_d2(r1,r2,r3,bb)
946 3072 : b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
947 768 : call dang_d1(r3,r2,r1,bb)
948 3136 : b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:)
949 : end do
950 :
951 : !third: dihedral values
952 2368 : do idihed=1,deloc%ndihed
953 2304 : i1 = deloc%dihedrals(1,1,idihed)
954 2304 : s1 = deloc%dihedrals(2,1,idihed)
955 : r1(:) = xcart(:,i1)+deloc%rshift(1,s1)*rprimd(:,1)&
956 : & +deloc%rshift(2,s1)*rprimd(:,2)&
957 9216 : & +deloc%rshift(3,s1)*rprimd(:,3)
958 2304 : i2 = deloc%dihedrals(1,2,idihed)
959 2304 : s2 = deloc%dihedrals(2,2,idihed)
960 : r2(:) = xcart(:,i2)+deloc%rshift(1,s2)*rprimd(:,1)&
961 : & +deloc%rshift(2,s2)*rprimd(:,2)&
962 9216 : & +deloc%rshift(3,s2)*rprimd(:,3)
963 2304 : i3 = deloc%dihedrals(1,3,idihed)
964 2304 : s3 = deloc%dihedrals(2,3,idihed)
965 : r3(:) = xcart(:,i3)+deloc%rshift(1,s3)*rprimd(:,1)&
966 : & +deloc%rshift(2,s3)*rprimd(:,2)&
967 9216 : & +deloc%rshift(3,s3)*rprimd(:,3)
968 2304 : i4 = deloc%dihedrals(1,4,idihed)
969 2304 : s4 = deloc%dihedrals(2,4,idihed)
970 : r4(:) = xcart(:,i4)+deloc%rshift(1,s4)*rprimd(:,1)&
971 : & +deloc%rshift(2,s4)*rprimd(:,2)&
972 9216 : & +deloc%rshift(3,s4)*rprimd(:,3)
973 : ! write(std_out,*) 'dihed ',idihed
974 : ! write(std_out,*) r1
975 : ! write(std_out,*) r2
976 : ! write(std_out,*) r3
977 : ! write(std_out,*) r4
978 :
979 2304 : iprim=iprim+1
980 2304 : call ddihedral_d1(r1,r2,r3,r4,bb)
981 9216 : b_matrix(iprim,3*(i1-1)+1:3*i1) = b_matrix(iprim,3*(i1-1)+1:3*i1) + bb(:)
982 2304 : call ddihedral_d2(r1,r2,r3,r4,bb)
983 9216 : b_matrix(iprim,3*(i2-1)+1:3*i2) = b_matrix(iprim,3*(i2-1)+1:3*i2) + bb(:)
984 2304 : call ddihedral_d2(r4,r3,r2,r1,bb)
985 9216 : b_matrix(iprim,3*(i3-1)+1:3*i3) = b_matrix(iprim,3*(i3-1)+1:3*i3) + bb(:)
986 2304 : call ddihedral_d1(r4,r3,r2,r1,bb)
987 9280 : b_matrix(iprim,3*(i4-1)+1:3*i4) = b_matrix(iprim,3*(i4-1)+1:3*i4) + bb(:)
988 : end do
989 :
990 64 : do icart=1,deloc%ncart
991 0 : iprim=iprim+1
992 : b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) = &
993 64 : & b_matrix(iprim,3*(deloc%carts(2,icart)-1)+deloc%carts(1,icart)) + one
994 : end do
995 :
996 : !DEBUG
997 : ! write (200,*) 'calc_b_matrix : b_matrix = '
998 : ! do iprim=1,deloc%ninternal
999 : ! do i1=1, 3*natom
1000 : ! write (200,'(E16.6,2x)',ADVANCE='NO') b_matrix(iprim,i1)
1001 : ! end do
1002 : ! write (200,*)
1003 : ! end do
1004 : !ENDDEBUG
1005 :
1006 64 : end subroutine calc_b_matrix
1007 : !!***
1008 :
1009 : !!****f* ABINIT/dbond_length_d1
1010 : !! NAME
1011 : !! dbond_length_d1
1012 : !!
1013 : !! FUNCTION
1014 : !!
1015 : !! SOURCE
1016 : !!
1017 :
1018 512 : subroutine dbond_length_d1(r1,r2,bb)
1019 :
1020 : !Arguments ------------------------------------
1021 : !arrays
1022 : real(dp),intent(in) :: r1(3),r2(3)
1023 : real(dp),intent(out) :: bb(3)
1024 :
1025 : !Local variables ------------------------------
1026 : !arrays
1027 : real(dp) :: rpt(3)
1028 :
1029 : !************************************************************************
1030 2048 : rpt(:) = r1(:)-r2(:)
1031 2048 : bb(:) = rpt(:)/bond_length(r1,r2)
1032 :
1033 512 : end subroutine dbond_length_d1
1034 : !!***
1035 :
1036 :
1037 : !!****f* ABINIT/dang_d1
1038 : !! NAME
1039 : !! dang_d1
1040 : !!
1041 : !! FUNCTION
1042 : !!
1043 : !! SOURCE
1044 : !!
1045 :
1046 1536 : subroutine dang_d1(r1,r2,r3,bb)
1047 :
1048 : !Arguments ------------------------------------
1049 : !arrays
1050 : real(dp),intent(in) :: r1(3),r2(3),r3(3)
1051 : real(dp),intent(out) :: bb(3)
1052 :
1053 : !Local variables ------------------------------
1054 : !scalars
1055 : real(dp) :: cos_ang,n1,n1232,n2,tmp
1056 : !arrays
1057 : real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3)
1058 :
1059 : !************************************************************************
1060 1536 : n1=bond_length(r1,r2)
1061 1536 : n2=bond_length(r3,r2)
1062 :
1063 6144 : rpt12(:) = r1(:)-r2(:)
1064 6144 : rpt32(:) = r3(:)-r2(:)
1065 :
1066 6144 : cos_ang = dot_product(rpt12,rpt32)/n1/n2
1067 1536 : if (cos_ang > one - epsilon(one)*two) then
1068 : cos_ang = one
1069 1536 : else if(cos_ang < -one + epsilon(one)*two) then
1070 0 : cos_ang = -one
1071 : end if
1072 :
1073 6144 : rpt(:) = rpt32(:)/n1/n2 - rpt12(:)*cos_ang/n1/n1
1074 :
1075 1536 : tmp = sqrt(one-cos_ang**2)
1076 1536 : bb(:) = zero
1077 1536 : if (tmp > epsilon(one)) then
1078 6144 : bb(:) = rpt(:) * (-one)/tmp
1079 : end if
1080 :
1081 : !TEST: version from MOLECULAR VIBRATIONS EB Wilson
1082 1536 : call acrossb(rpt12,rpt32,cp1232)
1083 6144 : n1232 = norm2(cp1232)
1084 6144 : rpt(:) = (cos_ang*rpt12(:)*n2/n1 - rpt32(:))/n1232
1085 1536 : if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then
1086 0 : write(std_out,*) 'Compare bb ang 1 : '
1087 0 : write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1088 : end if
1089 1536 : bb(:) = rpt(:)
1090 :
1091 1536 : end subroutine dang_d1
1092 : !!***
1093 :
1094 :
1095 : !!****f* ABINIT/dang_d2
1096 : !! NAME
1097 : !! dang_d2
1098 : !!
1099 : !! FUNCTION
1100 : !!
1101 : !! SOURCE
1102 : !!
1103 :
1104 768 : subroutine dang_d2(r1,r2,r3,bb)
1105 :
1106 : !Arguments ------------------------------------
1107 : !arrays
1108 : real(dp),intent(in) :: r1(3),r2(3),r3(3)
1109 : real(dp),intent(out) :: bb(3)
1110 :
1111 : !Local variables ------------------------------
1112 : !scalars
1113 : real(dp) :: cos_ang,n1,n1232,n2,tmp
1114 : !arrays
1115 : real(dp) :: cp1232(3),rpt(3),rpt12(3),rpt32(3)
1116 :
1117 : !************************************************************************
1118 768 : n1=bond_length(r1,r2)
1119 768 : n2=bond_length(r3,r2)
1120 :
1121 3072 : rpt12(:) = r1(:)-r2(:)
1122 3072 : rpt32(:) = r3(:)-r2(:)
1123 :
1124 3072 : cos_ang = dot_product(rpt12,rpt32)/n1/n2
1125 768 : if (cos_ang > one - epsilon(one)*two) then
1126 : cos_ang = one
1127 768 : else if(cos_ang < -one + epsilon(one)*two) then
1128 0 : cos_ang = -one
1129 : end if
1130 :
1131 : rpt(:) = -rpt32(:)/n1/n2 - rpt12(:)/n1/n2 &
1132 3072 : & + rpt12(:)*cos_ang/n1/n1 + rpt32(:)*cos_ang/n2/n2
1133 :
1134 768 : tmp = sqrt(one-cos_ang**2)
1135 768 : bb(:) = zero
1136 768 : if (tmp > tol12) then
1137 3072 : bb(:) = rpt(:) * (-one)/tmp
1138 : end if
1139 :
1140 : !TEST: version from MOLECULAR VIBRATIONS EB Wilson
1141 768 : call acrossb(rpt12,rpt32,cp1232)
1142 3072 : n1232 = norm2(cp1232)
1143 3072 : rpt(:) = ((n1-n2*cos_ang)*rpt12(:)/n1 + (n2-n1*cos_ang)*rpt32(:)/n2) / n1232
1144 768 : if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then
1145 0 : write(std_out,*) 'Compare bb ang 2 : '
1146 0 : write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1147 : end if
1148 768 : bb(:) = rpt(:)
1149 :
1150 768 : end subroutine dang_d2
1151 : !!***
1152 :
1153 : !!****f* ABINIT/ddihedral_d1
1154 : !! NAME
1155 : !! ddihedral_d1
1156 : !!
1157 : !! FUNCTION
1158 : !!
1159 : !! SOURCE
1160 : !!
1161 :
1162 4608 : subroutine ddihedral_d1(r1,r2,r3,r4,bb)
1163 :
1164 : !Arguments ------------------------------------
1165 : !arrays
1166 : real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
1167 : real(dp),intent(out) :: bb(3)
1168 :
1169 : !Local variables ------------------------------------
1170 : !scalars
1171 : real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,tmp
1172 : !arrays
1173 : real(dp) :: cp1232(3),cp32_1232(3),cp32_3432(3),cp3432(3),cpcp(3),rpt(3)
1174 : real(dp) :: rpt12(3),rpt32(3),rpt34(3)
1175 :
1176 : !******************************************************************
1177 18432 : rpt12(:) = r1(:)-r2(:)
1178 18432 : rpt32(:) = r3(:)-r2(:)
1179 18432 : rpt34(:) = r3(:)-r4(:)
1180 :
1181 4608 : call acrossb(rpt12,rpt32,cp1232)
1182 4608 : call acrossb(rpt34,rpt32,cp3432)
1183 :
1184 : !DEBUG
1185 : !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
1186 : !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
1187 : !ENDDEBUG
1188 :
1189 18432 : n1 = norm2(cp1232)
1190 18432 : n2 = norm2(cp3432)
1191 :
1192 18432 : cos_dihedral = dot_product(cp1232,cp3432)/n1/n2
1193 4608 : if (cos_dihedral > one - epsilon(one)*two) then
1194 : cos_dihedral = one
1195 3112 : else if(cos_dihedral < -one + epsilon(one)*two) then
1196 0 : cos_dihedral = -one
1197 : end if
1198 : !we use complementary of standard angle, so
1199 : !cos_dihedral = -cos_dihedral
1200 :
1201 4608 : call acrossb(cp1232,cp3432,cpcp)
1202 18432 : cpcp(:) = cpcp(:)/n1/n2
1203 : !we use complementary of standard angle, but sin is invariant
1204 32256 : sin_dihedral = -dot_product(cpcp,rpt32)/norm2(rpt32)
1205 4608 : dih_sign = one
1206 4608 : if (sin_dihedral < -epsilon(one)) then
1207 1556 : dih_sign = -one
1208 : end if
1209 :
1210 : !DEBUG
1211 : !write(std_out,'(a,3E16.6)') 'ddihedral_d1 : cos abs(sin) dih_sign= ',&
1212 : !& cos_dihedral,sin_dihedral,dih_sign
1213 : !ENDDEBUG
1214 :
1215 : !ddihedral_d1 = dih_sign* acos(cos_dihedral)
1216 4608 : call acrossb(rpt32,cp1232,cp32_1232)
1217 4608 : call acrossb(rpt32,cp3432,cp32_3432)
1218 :
1219 18432 : rpt(:) = cp32_3432(:)/n1/n2 - cp32_1232(:)/n1/n1 * cos_dihedral
1220 4608 : bb(:) = zero
1221 :
1222 : !DEBUG
1223 : !write(std_out,*) 'ddihedral_d1 cp1232 cp3432 = ',cp1232,cp3432,rpt32
1224 : !write(std_out,*) 'ddihedral_d1 cp32_1232 cp32_3432 = ',cp32_1232,cp32_3432,cos_dihedral,n1,n2
1225 : !write(std_out,*) 'ddihedral_d1 rpt = ',rpt
1226 : !ENDDEBUG
1227 :
1228 4608 : tmp = sqrt(one-cos_dihedral**2)
1229 4608 : if (tmp > tol12) then
1230 : ! we use complementary of standard angle, so cosine in acos has - sign,
1231 : ! and it appears for the derivative
1232 12448 : bb(:) = -dih_sign * rpt(:) * (-one) / tmp
1233 : else
1234 10472 : bb(:) = dih_sign * cp32_3432(:) / n1 / n2 / norm2(cp32_3432)
1235 : end if
1236 :
1237 : !TEST: version from MOLECULAR VIBRATIONS EB Wilson
1238 :
1239 18432 : n23 = norm2(rpt32)
1240 18432 : rpt(:) = cp1232(:)*n23/n1/n1
1241 : !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then
1242 : !write(std_out,*) 'Compare bb1 : '
1243 : !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1244 : !end if
1245 4608 : bb(:) = rpt(:)
1246 :
1247 4608 : end subroutine ddihedral_d1
1248 : !!***
1249 :
1250 : !!****f* ABINIT/ddihedral_d2
1251 : !! NAME
1252 : !! ddihedral_d2
1253 : !!
1254 : !! FUNCTION
1255 : !!
1256 : !! SOURCE
1257 : !!
1258 :
1259 4608 : subroutine ddihedral_d2(r1,r2,r3,r4,bb)
1260 :
1261 : !Arguments ------------------------------------
1262 : !arrays
1263 : real(dp),intent(in) :: r1(3),r2(3),r3(3),r4(3)
1264 : real(dp),intent(out) :: bb(3)
1265 :
1266 : !Local variables
1267 : !scalars
1268 : real(dp) :: cos_dihedral,dih_sign,n1,n2,n23,sin_dihedral,sp1232,sp3432,tmp
1269 : !arrays
1270 : real(dp) :: cp1232(3),cp1232_12(3),cp1232_34(3),cp32_1232(3),cp32_3432(3)
1271 : real(dp) :: cp3432(3),cp3432_12(3),cp3432_34(3),cpcp(3),rpt(3),rpt12(3)
1272 : real(dp) :: rpt32(3),rpt34(3)
1273 :
1274 : ! *************************************************************************
1275 18432 : rpt12(:) = r1(:)-r2(:)
1276 18432 : rpt32(:) = r3(:)-r2(:)
1277 18432 : rpt34(:) = r3(:)-r4(:)
1278 :
1279 4608 : call acrossb(rpt12,rpt32,cp1232)
1280 4608 : call acrossb(rpt34,rpt32,cp3432)
1281 :
1282 : !DEBUG
1283 : !write(std_out,*) ' cos_dihedral : cp1232 = ', cp1232
1284 : !write(std_out,*) ' cos_dihedral : cp3432 = ', cp3432
1285 : !ENDDEBUG
1286 :
1287 18432 : n1 = norm2(cp1232)
1288 18432 : n2 = norm2(cp3432)
1289 :
1290 18432 : cos_dihedral = dot_product(cp1232,cp3432)/n1/n2
1291 4608 : if (cos_dihedral > one - epsilon(one)*two) then
1292 : cos_dihedral = one
1293 3112 : else if(cos_dihedral < -one + epsilon(one)*two) then
1294 0 : cos_dihedral = -one
1295 : end if
1296 : !we use complementary of standard angle, so
1297 : !cos_dihedral = -cos_dihedral
1298 :
1299 4608 : call acrossb(cp1232,cp3432,cpcp)
1300 18432 : cpcp(:) = cpcp(:)/n1/n2
1301 : !we use complementary of standard angle, but sin is invariant
1302 32256 : sin_dihedral = -dot_product(cpcp,rpt32)/norm2(rpt32)
1303 4608 : dih_sign = one
1304 4608 : if (sin_dihedral < -tol12) then
1305 1536 : dih_sign = -one
1306 : end if
1307 :
1308 : !DEBUG
1309 : !write(std_out,'(a,3E16.6)') 'ddihedral_d2 : cos abs(sin) dih_sign= ',&
1310 : !& cos_dihedral,sin_dihedral,dih_sign
1311 : !ENDDEBUG
1312 :
1313 : !ddihedral_d2 = dih_sign* acos(cos_dihedral)
1314 4608 : call acrossb(rpt32,cp3432,cp32_3432)
1315 4608 : call acrossb(cp3432,rpt12,cp3432_12)
1316 4608 : call acrossb(cp1232,rpt34,cp1232_34)
1317 :
1318 4608 : call acrossb(rpt32,cp1232,cp32_1232)
1319 4608 : call acrossb(cp1232,rpt12,cp1232_12)
1320 4608 : call acrossb(cp3432,rpt34,cp3432_34)
1321 :
1322 : rpt(:) = -(cp32_3432(:) + cp3432_12(:) + cp1232_34(:))/n1/n2 &
1323 18432 : & +cos_dihedral*(cp32_1232(:)/n1/n1 + cp1232_12(:)/n1/n1 + cp3432_34(:)/n2/n2)
1324 4608 : bb(:) = zero
1325 4608 : tmp = sqrt(one-cos_dihedral**2)
1326 4608 : if (tmp > tol12) then
1327 : ! we use complementary of standard angle, so cosine in acos has - sign,
1328 : ! and it appears for derivative
1329 12448 : bb(:) = -dih_sign * rpt(:) * (-one) / tmp
1330 : else
1331 : bb(:) = dih_sign * cos_dihedral * &
1332 : & ( cp32_1232(:)/n1/n1/sqrt(cp32_1232(1)**2+cp32_1232(2)**2+cp32_1232(3)**2) &
1333 : & +cp1232_12(:)/n1/n1/sqrt(cp1232_12(1)**2+cp1232_12(2)**2+cp1232_12(3)**2) &
1334 5984 : & +cp3432_34(:)/n2/n2/sqrt(cp3432_34(1)**2+cp3432_34(2)**2+cp3432_34(3)**2) )
1335 : end if
1336 :
1337 : !TEST: version from MOLECULAR VIBRATIONS EB Wilson p. 61
1338 18432 : n23 = norm2(rpt32)
1339 18432 : sp1232 = dot_product(rpt12,rpt32)
1340 18432 : sp3432 = dot_product(rpt34,rpt32)
1341 :
1342 18432 : rpt(:) = -cp1232(:)*(n23-sp1232/n23)/n1/n1 - cp3432(:)*sp3432/n23/n2/n2
1343 : !if (abs(bb(1)-rpt(1))+abs(bb(2)-rpt(2))+abs(bb(3)-rpt(3)) > tol10) then
1344 : !write(std_out,*) 'Compare bb2 : '
1345 : !write(std_out,*) bb(:), rpt(:), bb(:)-rpt(:)
1346 : !write(std_out,*) -cp1232(:)*(n23-sp1232/n23)/n1/n1, -cp3432(:)*sp3432/n23/n2/n2
1347 : !end if
1348 4608 : bb(:) = rpt(:)
1349 :
1350 4608 : end subroutine ddihedral_d2
1351 : !!***
1352 :
1353 : !!****f* ABINIT/xcart2deloc
1354 : !! NAME
1355 : !! xcart2deloc
1356 : !!
1357 : !! FUNCTION
1358 : !! Calculate values of delocalized coordinates as a function of
1359 : !! cartesian ones. First primitive internals, then B matrix,
1360 : !! then F, then U then delocalized internals.
1361 : !!
1362 : !! INPUTS
1363 : !! deloc <type(delocint)>=Important variables for pred_delocint
1364 : !! |
1365 : !! | nang = Number of angles
1366 : !! | nbond = Number of bonds
1367 : !! | ncart = Number of cartesian directions
1368 : !! | (used for constraints)
1369 : !! | ndihed = Number of dihedrals
1370 : !! | nrshift = Dimension of rshift
1371 : !! | ninternal= Number of internal coordinates
1372 : !! | ninternal=nbond+nang+ndihed+ncart
1373 : !! |
1374 : !! | angs(2,3,nang) = Indexes to characterize angles
1375 : !! | bonds(2,2,nbond)= For a bond between iatom and jatom
1376 : !! | bonds(1,1,nbond) = iatom
1377 : !! | bonds(2,1,nbond) = icenter
1378 : !! | bonds(1,2,nbond) = jatom
1379 : !! | bonds(2,2,nbond) = irshift
1380 : !! | carts(2,ncart) = Index of total primitive internal,
1381 : !! | and atom (carts(2,:))
1382 : !! | dihedrals(2,4,ndihed)= Indexes to characterize dihedrals
1383 : !! |
1384 : !! | rshift(3,nrshift)= Shift in xred that must be done to find
1385 : !! | all neighbors of a given atom within a
1386 : !! | given number of neighboring shells
1387 : !! natom = Number of atoms
1388 : !! rprimd(3,3) = Dimensional real space primitive translations
1389 : !! (bohr)
1390 : !! xcart(3,natom) = Cartesian coordinates of atoms (bohr)
1391 : !!
1392 : !! OUTPUT
1393 : !! bt_inv_matrix(3*(natom-1),3*natom) = Inverse of B^{T} matrix
1394 : !! deloc_int(3*(natom-1)) = Delocalized internal coordinates
1395 : !! prim_int(ninternal) = Primitive internal coordinates
1396 : !!
1397 : !! SIDE EFFECTS
1398 : !! u_matrix(ninternal,3*(natom-1)) = Eigenvectors of BB^T matrix
1399 : !!
1400 : !! NOTES
1401 : !!
1402 : !! SOURCE
1403 :
1404 64 : subroutine xcart2deloc(deloc,natom,rprimd,xcart,bt_inv_matrix,u_matrix,deloc_int,prim_int)
1405 :
1406 : !Arguments ------------------------------------
1407 : !scalars
1408 : integer,intent(in) :: natom
1409 : type(delocint),intent(in) :: deloc
1410 : !arrays
1411 : real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
1412 : real(dp),intent(inout) :: u_matrix(deloc%ninternal,3*(natom-1))
1413 : real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom)
1414 : real(dp),intent(out) :: deloc_int(3*(natom-1))
1415 : real(dp),intent(out) :: prim_int(deloc%ninternal)
1416 :
1417 : !Local variables-------------------------------
1418 : !scalars
1419 : integer :: ii
1420 : logical :: DEBUG=.FALSE.
1421 : !arrays
1422 128 : real(dp) :: b_matrix(deloc%ninternal,3*natom)
1423 :
1424 : ! ******************************************************************
1425 :
1426 64 : call calc_prim_int(deloc,natom,rprimd,xcart,prim_int)
1427 64 : if (DEBUG)then
1428 0 : write(std_out,*) 'Primitive Internals'
1429 0 : do ii=1,deloc%ninternal
1430 0 : write(std_out,*) prim_int(ii)
1431 : end do
1432 : end if
1433 :
1434 64 : call calc_b_matrix(deloc,natom,rprimd,xcart,b_matrix)
1435 64 : if (DEBUG)then
1436 0 : write(std_out,*) 'B Matrix'
1437 0 : do ii=1,deloc%ninternal
1438 0 : write(std_out,*) b_matrix(:,ii)
1439 : end do
1440 : end if
1441 :
1442 : call calc_btinv_matrix(b_matrix,natom,deloc%ninternal,&
1443 64 : & bt_inv_matrix,u_matrix)
1444 64 : if (DEBUG)then
1445 0 : write(std_out,*) 'BT Inverse Matrix'
1446 0 : do ii=1,3*natom
1447 0 : write(std_out,*) bt_inv_matrix(:,ii)
1448 : end do
1449 : end if
1450 :
1451 : !calculate value of delocalized internals
1452 :
1453 : call dgemv('T',deloc%ninternal,3*(natom-1),one,&
1454 64 : & u_matrix,deloc%ninternal,prim_int,1,zero,deloc_int,1)
1455 :
1456 64 : end subroutine xcart2deloc
1457 : !!***
1458 :
1459 :
1460 : !!****f* ABINIT/calc_btinv_matrix
1461 : !! NAME
1462 : !! calc_btinv_matrix
1463 : !!
1464 : !! FUNCTION
1465 : !!
1466 : !! INPUTS
1467 : !!
1468 : !! OUTPUT
1469 : !!
1470 : !! NOTES
1471 : !! bt_inv_matrix is inverse transpose of the delocalized
1472 : !! coordinate B matrix. b_matrix is the primitive internal B matrix
1473 : !!
1474 : !! SOURCE
1475 :
1476 64 : subroutine calc_btinv_matrix(b_matrix,natom,ninternal,bt_inv_matrix,u_matrix)
1477 :
1478 : !Arguments ------------------------------------
1479 : integer,intent(in) :: ninternal,natom
1480 : real(dp),intent(in) :: b_matrix(ninternal,3*natom)
1481 : real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom)
1482 : real(dp),intent(inout) :: u_matrix(ninternal,3*(natom-1))
1483 :
1484 : !Local variables ------------------------------------
1485 : !scalars
1486 : integer :: ii,info,lwork
1487 : !arrays
1488 128 : real(dp) :: f_eigs(3*natom),f_matrix(3*natom,3*natom)
1489 128 : real(dp) :: s_matrix(3*natom,3*natom)
1490 128 : real(dp) :: s_red(3*natom,3*(natom-1))
1491 128 : real(dp) :: u_matrix_old(ninternal,3*(natom-1))
1492 64 : real(dp),allocatable :: work(:)
1493 :
1494 : !******************************************************************
1495 :
1496 : !f matrix = B^{T} B
1497 : call dgemm('T','N',3*natom,3*natom,ninternal,one,&
1498 64 : & b_matrix,ninternal,b_matrix,ninternal,zero,f_matrix,3*natom)
1499 :
1500 64 : lwork = max(1,3*3*natom-1)
1501 192 : ABI_MALLOC(work,(lwork))
1502 2752 : s_matrix(:,:) = f_matrix(:,:)
1503 :
1504 64 : call dsyev('V','L',3*natom,s_matrix,3*natom,f_eigs,work,lwork,info)
1505 :
1506 64 : ABI_FREE(work)
1507 :
1508 64 : if (abs(f_eigs(1)) + abs(f_eigs(2)) + abs(f_eigs(3)) > tol10 ) then
1509 0 : write(std_out,*) 'Error: 3 lowest eigenvalues are not zero'
1510 0 : write(std_out,*) ' internal coordinates do NOT span the full degrees of freedom !'
1511 0 : write(std_out,'(6E16.6)') f_eigs
1512 0 : ABI_ERROR("Aborting now")
1513 : end if
1514 64 : if ( abs(f_eigs(4)) < tol10 ) then
1515 0 : write(std_out,*) 'Error: fourth eigenvalue is zero'
1516 0 : write(std_out,*) ' internal coordinates do NOT span the full degrees of freedom !'
1517 0 : write(std_out,'(6E16.6)') f_eigs
1518 0 : ABI_ERROR("Aborting now")
1519 : end if
1520 :
1521 : !calculate U matrix from U = B * S_red * lambda^{-1/2}
1522 256 : do ii=1,3*(natom-1)
1523 1408 : s_red(:,ii) = s_matrix(:,ii+3)/sqrt(f_eigs(ii+3))
1524 : end do
1525 :
1526 10240 : u_matrix_old(:,:) = u_matrix(:,:)
1527 :
1528 : call dgemm('N','N',ninternal,3*(natom-1),3*natom,one,&
1529 64 : & b_matrix,ninternal,s_red,3*natom,zero,u_matrix,ninternal)
1530 :
1531 :
1532 : !align eigenvectors, to preserve a form of continuity in convergences
1533 : !!!! eigenvalues are no longer in increasing order!!! but only s_red is reordered
1534 : !so that btinv is correct.
1535 64 : call align_u_matrices(natom,ninternal,u_matrix,u_matrix_old,s_matrix,f_eigs)
1536 :
1537 : !calculate B_deloc^{-1} matrix for transformation of forces to deloc coord.
1538 : !(B^{T}_deloc)^{-1} = (B_deloc B^{T}_deloc)^{-1} B_deloc = lambda^{-3/2} S^{T} F
1539 : != ( S lambda^{3/2} )^{T} F
1540 :
1541 : !! DEFINITION
1542 : !! real(dp),intent(out) :: bt_inv_matrix(3*(natom-1),3*natom)
1543 :
1544 : !even better: B_deloc^{-1} = lambda^{-1/2} S^{T}
1545 256 : do ii=1,3*(natom-1)
1546 : ! s_red(:,ii) = s_matrix(:,ii+3)*sqrt(f_eigs(ii+3))
1547 1408 : bt_inv_matrix(ii,:) = s_matrix(:,ii+3)/sqrt(f_eigs(ii+3))
1548 : end do
1549 :
1550 64 : end subroutine calc_btinv_matrix
1551 : !!***
1552 :
1553 : !!****f* ABINIT/align_u_matrices
1554 : !! NAME
1555 : !! align_u_matrices
1556 : !!
1557 : !! FUNCTION
1558 : !!
1559 : !! INPUTS
1560 : !!
1561 : !! OUTPUT
1562 : !!
1563 : !! SOURCE
1564 :
1565 64 : subroutine align_u_matrices(natom,ninternal,u_matrix,u_matrix_old,s_matrix,f_eigs)
1566 :
1567 : !Arguments ------------------------------------
1568 : !scalars
1569 : integer,intent(in) :: ninternal,natom
1570 : !arrays
1571 : real(dp),intent(in) :: u_matrix_old(ninternal,3*(natom-1))
1572 : real(dp),intent(inout) :: f_eigs(3*natom)
1573 : real(dp),intent(inout) :: s_matrix(3*natom,3*natom)
1574 : real(dp),intent(inout) :: u_matrix(ninternal,3*(natom-1))
1575 :
1576 : !Local variables ------------------------------
1577 : !scalars
1578 : integer :: ii,iint1,imax
1579 : real(dp) :: ss
1580 : !arrays
1581 128 : integer :: eigv_flag(3*(natom-1)),eigv_ind(3*(natom-1))
1582 128 : real(dp) :: tmps(3*natom,3*natom)
1583 128 : real(dp) :: tmpu(ninternal,3*(natom-1))
1584 64 : real(dp) :: tmpf(3*natom)
1585 :
1586 : !******************************************************************
1587 :
1588 256 : eigv_flag(:) = 0
1589 256 : eigv_ind(:) = 0
1590 :
1591 : !just permit a change in sign
1592 256 : do iint1=1,3*(natom-1)
1593 : ss = zero
1594 10176 : do ii=1,ninternal
1595 10176 : ss = ss + u_matrix_old(ii,iint1)*u_matrix(ii,iint1)
1596 : end do
1597 192 : if (ss < -tol12) then
1598 0 : imax = -iint1
1599 : else
1600 : imax = iint1
1601 : end if
1602 192 : eigv_ind(iint1) = imax
1603 256 : eigv_flag(abs(imax)) = 1
1604 : end do
1605 :
1606 10240 : tmpu(:,:) = u_matrix
1607 2752 : tmps(:,:) = s_matrix
1608 448 : tmpf(:) = f_eigs
1609 : !exchange eigenvectors...
1610 256 : do iint1=1,3*(natom-1)
1611 192 : ss = one
1612 192 : if (eigv_ind(iint1) < 0) ss = -one
1613 :
1614 192 : imax = abs(eigv_ind(iint1))
1615 :
1616 10176 : tmpu(:,imax) = ss*u_matrix(:,iint1)
1617 :
1618 1344 : tmps(:,imax+3) = ss*s_matrix(:,iint1+3)
1619 :
1620 256 : tmpf(imax+3) = f_eigs(iint1+3)
1621 : end do
1622 :
1623 10240 : u_matrix(:,:) = tmpu(:,:)
1624 2752 : s_matrix(:,:) = tmps(:,:)
1625 448 : f_eigs(:) = tmpf(:)
1626 :
1627 64 : end subroutine align_u_matrices
1628 : !!***
1629 :
1630 : !!****f* ABINIT/xfh_recover_deloc
1631 : !! NAME
1632 : !! xfh_recover_deloc
1633 : !!
1634 : !! FUNCTION
1635 : !! Update the contents of the history xfhist taking values
1636 : !! from xred, acell, rprim, gred_corrected and strten
1637 : !!
1638 : !! INPUTS
1639 : !!
1640 : !! OUTPUT
1641 : !!
1642 : !! SOURCE
1643 :
1644 0 : subroutine xfh_recover_deloc(ab_xfh,ab_mover,acell,cycl_main,&
1645 0 : & gred,hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,&
1646 0 : & vout,vout_prev,xred,deloc,deloc_int,deloc_gred,btinv,gprimd,prim_int,&
1647 0 : & u_matrix)
1648 :
1649 : !Arguments ------------------------------------
1650 : !scalars
1651 :
1652 : integer,intent(in) :: ndim
1653 : integer,intent(out) :: cycl_main
1654 : real(dp),intent(inout) :: ucvol,ucvol0
1655 : type(ab_xfh_type),intent(inout) :: ab_xfh
1656 : type(abimover),intent(in) :: ab_mover
1657 : ! DELOCINT specials
1658 : type(delocint),intent(in) :: deloc
1659 :
1660 : !arrays
1661 : real(dp),intent(inout) :: acell(3)
1662 : real(dp),intent(inout) :: hessin(:,:)
1663 : real(dp),intent(inout) :: xred(3,ab_mover%natom)
1664 : real(dp),intent(inout) :: rprim(3,3)
1665 : real(dp),intent(inout) :: rprimd0(3,3)
1666 : real(dp),intent(inout) :: gred(3,ab_mover%natom)
1667 : real(dp),intent(inout) :: strten(6)
1668 : real(dp),intent(inout) :: vin(:)
1669 : real(dp),intent(inout) :: vin_prev(:)
1670 : real(dp),intent(inout) :: vout(:)
1671 : real(dp),intent(inout) :: vout_prev(:)
1672 : ! DELOCINT specials
1673 : real(dp),intent(inout) :: deloc_gred(3*(ab_mover%natom-1))
1674 : real(dp),intent(inout) :: deloc_int(3*(ab_mover%natom-1))
1675 : real(dp),intent(inout) :: btinv(3*(ab_mover%natom-1),3*ab_mover%natom)
1676 : real(dp),intent(inout) :: prim_int(:),u_matrix(:,:),gprimd(3,3)
1677 :
1678 : !Local variables-------------------------------
1679 : !scalars
1680 : integer :: ixfh
1681 0 : real(dp) :: xcart(3,ab_mover%natom)
1682 :
1683 : !*********************************************************************
1684 :
1685 0 : if(ab_xfh%nxfh/=0)then
1686 : ! Loop over previous time steps
1687 0 : do ixfh=1,ab_xfh%nxfh
1688 :
1689 : ! For that time step, get new (x,f) from xfhist
1690 0 : xred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom ,1,ixfh)
1691 0 : rprim(1:3,1:3)=ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh)
1692 0 : acell(:) =ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh)
1693 0 : gred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh)
1694 : ! This use of results_gs is unusual
1695 0 : strten(1:3) =ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh)
1696 0 : strten(4:6) =ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh)
1697 :
1698 : ! !DEBUG
1699 : ! write (ab_out,*) '---READ FROM XFHIST---'
1700 :
1701 : ! write (ab_out,*) 'XRED'
1702 : ! do kk=1,ab_mover%natom
1703 : ! write (ab_out,*) xred(:,kk)
1704 : ! end do
1705 : ! write (ab_out,*) 'FRED'
1706 : ! do kk=1,ab_mover%natom
1707 : ! write (ab_out,*) gred(:,kk)
1708 : ! end do
1709 : ! write(ab_out,*) 'RPRIM'
1710 : ! do kk=1,3
1711 : ! write(ab_out,*) rprim(:,kk)
1712 : ! end do
1713 : ! write(ab_out,*) 'ACELL'
1714 : ! write(ab_out,*) acell(:)
1715 : ! !DEBUG
1716 :
1717 : ! Convert input xred (reduced coordinates) to xcart (cartesian)
1718 0 : call xred2xcart(ab_mover%natom,rprimd0,xcart,xred)
1719 : ! Convert input coordinates in Delocalized internals
1720 : call xcart2deloc(deloc,ab_mover%natom,rprimd0,xcart,&
1721 0 : & btinv,u_matrix,deloc_int,prim_int)
1722 : ! Convert forces to delocalized coordinates for next step
1723 0 : call gred2gdeloc(btinv,deloc_gred,gred,ab_mover%natom,gprimd)
1724 :
1725 : ! Transfer it in vin, vout
1726 : call xfpack_x2vin(acell,ab_mover%natom-1,&
1727 : & ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,&
1728 0 : & ab_mover%symrel,ucvol,ucvol0,vin,deloc_int)
1729 : call xfpack_f2vout(deloc_gred,ab_mover%natom-1,&
1730 : & ndim,ab_mover%optcell,ab_mover%strtarget,strten,&
1731 0 : & ucvol,vout)
1732 : ! Get old time step, if any, and update inverse hessian
1733 0 : if(ixfh/=1)then
1734 0 : xred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom,1,ixfh-1)
1735 : rprim(1:3,1:3)=&
1736 0 : & ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh-1)
1737 0 : acell(:)=ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh-1)
1738 0 : gred(:,:)=ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh-1)
1739 : ! This use of results_gs is unusual
1740 0 : strten(1:3)=ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh-1)
1741 0 : strten(4:6)=ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh-1)
1742 :
1743 : ! Convert input xred (reduced coordinates) to xcart (cartesian)
1744 0 : call xred2xcart(ab_mover%natom,rprimd0,xcart,xred)
1745 : ! Convert input coordinates in Delocalized internals
1746 : call xcart2deloc(deloc,ab_mover%natom,rprimd0,xcart,&
1747 0 : & btinv,u_matrix,deloc_int,prim_int)
1748 : ! Convert forces to delocalized coordinates for next step
1749 0 : call gred2gdeloc(btinv,deloc_gred,gred,ab_mover%natom,gprimd)
1750 :
1751 : ! Tranfer it in vin_prev, vout_prev
1752 : call xfpack_x2vin(acell,ab_mover%natom-1,&
1753 : & ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,&
1754 0 : & ab_mover%symrel,ucvol,ucvol0,vin_prev,deloc_int)
1755 : call xfpack_f2vout(deloc_gred,ab_mover%natom-1,&
1756 : & ndim,ab_mover%optcell,ab_mover%strtarget,strten,&
1757 0 : & ucvol,vout_prev)
1758 :
1759 : ! write(ab_out,*) 'Hessian matrix before update',ndim,'x',ndim
1760 : ! write(ab_out,*) 'ixfh=',ixfh
1761 : ! do kk=1,ndim
1762 : ! do jj=1,ndim,3
1763 : ! if (jj+2<=ndim)then
1764 : ! write(ab_out,*) jj,hessin(jj:jj+2,kk)
1765 : ! else
1766 : ! write(ab_out,*) jj,hessin(jj:ndim,kk)
1767 : ! end if
1768 : ! end do
1769 : ! end do
1770 :
1771 : call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom-1,ndim,&
1772 0 : & vin,vin_prev,vout,vout_prev)
1773 :
1774 : ! !DEBUG
1775 : ! write(ab_out,*) 'Hessian matrix after update',ndim,'x',ndim
1776 : ! do kk=1,ndim
1777 : ! do jj=1,ndim,3
1778 : ! if (jj+2<=ndim)then
1779 : ! write(ab_out,*) jj,hessin(jj:jj+2,kk)
1780 : ! else
1781 : ! write(ab_out,*) jj,hessin(jj:ndim,kk)
1782 : ! end if
1783 : ! end do
1784 : ! end do
1785 : ! !DEBUG
1786 :
1787 : end if !if(ab_xfh%nxfh/=0)
1788 :
1789 : ! End loop over previous time steps
1790 : end do
1791 :
1792 : ! The hessian has been generated,
1793 : ! as well as the latest vin and vout
1794 : ! so will cycle the main loop
1795 0 : cycl_main=1
1796 :
1797 : end if
1798 :
1799 0 : end subroutine xfh_recover_deloc
1800 : !!***
1801 :
1802 1 : end module m_pred_delocint
1803 : !!***
|
