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