Line data Source code
1 : !!****m* ABINIT/m_xfpack
2 : !! NAME
3 : !! m_xfpack
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (XG, MJV, DCA, 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_xfpack
22 :
23 : use defs_basis
24 : use m_errors
25 : use m_abicore
26 : use m_abimover
27 :
28 : use m_matrix, only : matr3inv
29 : use m_geometry, only : mkradim, mkrdim, metric, strainsym
30 : use m_results_gs , only : results_gs_type
31 : use m_bfgs, only : hessupdt
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: xfpack_vin2x
39 : public :: xfpack_x2vin
40 : public :: xfpack_f2vout
41 : public :: xfh_recover_new
42 : public :: xfh_update
43 : !!***
44 :
45 : contains
46 : !!***
47 :
48 : !!****f* ABINIT/xfpack_vin2x
49 : !! NAME
50 : !! xfpack_vin2x
51 : !!
52 : !! FUNCTION
53 : !! Old option=2, transfer vin to xred, acell and rprim
54 : !!
55 : !! INPUTS
56 : !! acell0(3)=reference length scales of primitive translations (bohr), needed for some values of optcell.
57 : !! natom=number of atoms in cell
58 : !! ndim=dimension of vin array
59 : !! nsym=order of group.
60 : !! rprimd0(3,3)=reference real space primitive translations,
61 : !! needed for some values of optcell.
62 : !! optcell=option for the optimisation of the unit cell. Described in abinit_help.
63 : !! Depending on its value, different part of acell and rprim
64 : !! are contained in vin.
65 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
66 : !! ucvol=unit cell volume (bohr^3), needed for some values of optcell.
67 : !! ucvol0=reference unit cell volume (bohr^3), needed for some values of optcell.
68 : !!
69 : !! OUTPUT
70 : !! (see side effects)
71 : !!
72 : !! SIDE EFFECTS
73 : !! Input/Output variables
74 : !! acell(3)=length scales of primitive translations (bohr)
75 : !! rprim(3,3)=dimensionless real space primitive translations
76 : !! vin(ndim)=vector that contains xred and some quantity derived
77 : !! from acell and rprim, depending on the value of optcell.
78 : !! xred(3,natom)=reduced dimensionless atomic coordinates
79 : !!
80 : !! SOURCE
81 :
82 1009 : subroutine xfpack_vin2x(acell,acell0,natom,ndim,nsym,optcell,&
83 1009 : & rprim,rprimd0,symrel,ucvol,ucvol0,vin,xred)
84 :
85 : !Arguments ------------------------------------
86 : !scalars
87 : integer,intent(in) :: natom,ndim,nsym,optcell
88 : real(dp),intent(in) :: ucvol0
89 : real(dp),intent(out) :: ucvol
90 : !arrays
91 : integer,intent(in) :: symrel(3,3,nsym)
92 : real(dp),intent(in) :: acell0(3),rprimd0(3,3)
93 : real(dp),intent(inout) :: acell(3),rprim(3,3)
94 : real(dp),intent(in) :: vin(ndim)
95 : real(dp),intent(out) :: xred(3,natom)
96 :
97 : !Local variables-------------------------------
98 : !scalars
99 : integer :: ii,jj,kk
100 : real(dp) :: scale
101 : character(len=500) :: message
102 : logical :: equal=.TRUE.
103 : !arrays
104 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3)
105 : real(dp) :: rprimd_symm(3,3),scaling(3,3)
106 :
107 : ! *************************************************************************
108 :
109 : !!DEBUG
110 : !write(ab_out,*) ''
111 : !write(ab_out,*) 'xfpack_vin2x'
112 : !write(ab_out,*) 'natom=',natom
113 : !write(ab_out,*) 'ndim=',ndim
114 : !write(ab_out,*) 'nsym=',nsym
115 : !write(ab_out,*) 'optcell=',optcell
116 : !write(ab_out,*) 'ucvol=',ucvol
117 : !write(ab_out,*) 'xred='
118 : !do kk=1,natom
119 : !write(ab_out,*) xred(:,kk)
120 : !end do
121 : !write(ab_out,*) 'VECTOR INPUT (vin) xfpack_vin2x INPUT'
122 : !do ii=1,ndim,3
123 : !if (ii+2<=ndim)then
124 : !write(ab_out,*) ii,vin(ii:ii+2)
125 : !else
126 : !write(ab_out,*) ii,vin(ii:ndim)
127 : !end if
128 : !end do
129 : !!DEBUG
130 :
131 :
132 : !##########################################################
133 : !### 1. Test for compatible ndim
134 :
135 1009 : if(optcell==0 .and. ndim/=3*natom)then
136 : write(message,'(a,a,a,i4,a,i4,a)' )&
137 0 : & ' When optcell=0, ndim MUST be equal to 3*natom,',ch10,&
138 0 : & ' while ndim=',ndim,' and 3*natom=',3*natom,'.'
139 0 : ABI_BUG(messagE)
140 : end if
141 :
142 1009 : if( (optcell==1) .and. ndim/=3*natom+1)then
143 : write(message,'(a,a,a,i4,a,i4,a)' )&
144 0 : & ' When optcell=1 ndim MUST be equal to 3*natom+1,',ch10,&
145 0 : & ' while ndim=',ndim,' and 3*natom+1=',3*natom+1,'.'
146 0 : ABI_BUG(message)
147 : end if
148 :
149 1009 : if( (optcell==2 .or. optcell==3) .and. ndim/=3*natom+6) then
150 : write(message,'(a,a,a,i4,a,i4,a)' )&
151 0 : & ' When optcell=2 or 3, ndim MUST be equal to 3*natom+6,',ch10,&
152 0 : & ' while ndim=',ndim,' and 3*natom+6=',3*natom+6,'.'
153 0 : ABI_BUG(message)
154 : end if
155 :
156 1009 : if( optcell>=4 .and. ndim/=3*natom+3)then
157 : write(message,'(a,a,a,i4,a,i4,a)' )&
158 0 : & ' When optcell=4,5,6,7,8 or 9, ndim MUST be equal to 3*natom+3,',ch10,&
159 0 : & ' while ndim=',ndim,' and 3*natom+3=',3*natom+3,'.'
160 0 : ABI_BUG(message)
161 : end if
162 :
163 : !##########################################################
164 : !### 3. option=2, transfer vin to xred, acell and rprim
165 :
166 : !Get xred, and eventually acell and rprim from vin
167 3027 : xred(:,:)=reshape( vin(1:3*natom), (/3,natom/) )
168 :
169 1009 : if(optcell==1)then
170 :
171 : ! acell(:)=acell0(:)*vin(3*natom+1)/(ucvol0**third)
172 632 : acell(:)=acell0(:)*vin(3*natom+1)
173 :
174 851 : else if (optcell>=2)then
175 :
176 160 : scaling(:,:)=0.0_dp
177 160 : scaling(1,1)=1.0_dp ; scaling(2,2)=1.0_dp ; scaling(3,3)=1.0_dp
178 :
179 : if(optcell==2 .or. optcell==3)then
180 129 : scaling(1,1)=vin(3*natom+1)
181 129 : scaling(2,2)=vin(3*natom+2)
182 129 : scaling(3,3)=vin(3*natom+3)
183 129 : scaling(2,3)=vin(3*natom+4) ; scaling(3,2)=vin(3*natom+4)
184 129 : scaling(1,3)=vin(3*natom+5) ; scaling(3,1)=vin(3*natom+5)
185 129 : scaling(1,2)=vin(3*natom+6) ; scaling(2,1)=vin(3*natom+6)
186 : else if(optcell==4)then
187 9 : scaling(1,1)=vin(3*natom+1)
188 9 : if (abs(scaling(1,1) - 1.0_dp) <1.E-14) scaling(1,1)=1.0_dp
189 9 : scaling(3,1)=vin(3*natom+2)
190 9 : if (abs(scaling(3,1)) <1.E-14) scaling(3,1)=0.0_dp
191 9 : scaling(2,1)=vin(3*natom+3)
192 9 : if (abs(scaling(2,1)) <1.E-14) scaling(2,1)=0.0_dp
193 : else if(optcell==5)then
194 0 : scaling(2,2)=vin(3*natom+1)
195 0 : if (abs(scaling(2,2) - 1.0_dp) <1.E-14) scaling(2,2)=1.0_dp
196 0 : scaling(3,2)=vin(3*natom+2)
197 0 : if (abs(scaling(3,2)) <1.E-14) scaling(3,2)=0.0_dp
198 0 : scaling(1,2)=vin(3*natom+3)
199 0 : if (abs(scaling(1,2)) <1.E-14) scaling(1,2)=0.0_dp
200 : else if(optcell==6)then
201 5 : scaling(3,3)=vin(3*natom+1)
202 5 : if (abs(scaling(3,3) - 1.0_dp) <1.E-14) scaling(3,3)=1.0_dp
203 5 : scaling(2,3)=vin(3*natom+2)
204 5 : if (abs(scaling(2,3)) <1.E-14) scaling(2,3)=0.0_dp
205 5 : scaling(1,3)=vin(3*natom+3)
206 5 : if (abs(scaling(1,3)) <1.E-14) scaling(1,3)=0.0_dp
207 : else if(optcell==7)then
208 13 : scaling(2,2)=vin(3*natom+2) ; scaling(3,3)=vin(3*natom+3)
209 13 : scaling(2,3)=vin(3*natom+1) ; scaling(3,2)=vin(3*natom+1)
210 : else if(optcell==8)then
211 0 : scaling(1,1)=vin(3*natom+1) ; scaling(3,3)=vin(3*natom+3)
212 0 : scaling(1,3)=vin(3*natom+2) ; scaling(3,1)=vin(3*natom+2)
213 : else if(optcell==9)then
214 4 : scaling(1,1)=vin(3*natom+1) ; scaling(2,2)=vin(3*natom+2)
215 4 : scaling(1,2)=vin(3*natom+3) ; scaling(2,1)=vin(3*natom+3)
216 : end if
217 160 : if(optcell<=3 .or. optcell>=7)then
218 584 : do ii=1,3
219 1898 : do jj=1,3
220 1314 : rprimd(ii,jj)=0.0_dp
221 5694 : do kk=1,3
222 5256 : rprimd(ii,jj)=rprimd(ii,jj)+scaling(ii,kk)*rprimd0(kk,jj)
223 : end do
224 : end do
225 : end do
226 : ! for optcell=4,5,6, implementing search for all 3 components of the vector to be relaxed according to Eq.10 of J. Chem. Phys.
227 : ! 136, 074103 (2012), i.e. search direction given by rprimd0 * stress
228 14 : else if(optcell==4)then
229 36 : rprimd(:,2) = rprimd0(:,2)
230 36 : rprimd(:,3) = rprimd0(:,3)
231 36 : rprimd(:,1) = 0.0_dp
232 36 : do ii=1,3
233 117 : do kk=1,3
234 108 : rprimd(ii,1) = rprimd(ii,1) + scaling(kk,1)*rprimd0(ii,kk)
235 : end do
236 : end do
237 5 : else if(optcell==5)then
238 0 : rprimd(:,1) = rprimd0(:,1)
239 0 : rprimd(:,3) = rprimd0(:,3)
240 0 : rprimd(:,2) = 0.0_dp
241 0 : do ii=1,3
242 0 : do kk=1,3
243 0 : rprimd(ii,2) = rprimd(ii,2) + scaling(kk,2)*rprimd0(ii,kk)
244 : end do
245 : end do
246 : else if(optcell==6)then
247 20 : rprimd(:,1) = rprimd0(:,1)
248 20 : rprimd(:,2) = rprimd0(:,2)
249 20 : rprimd(:,3) = 0.0_dp
250 20 : do ii=1,3
251 65 : do kk=1,3
252 60 : rprimd(ii,3) = rprimd(ii,3) + scaling(kk,3)*rprimd0(ii,kk)
253 : end do
254 : end do
255 : end if
256 :
257 : ! Rescale if the volume must be preserved
258 160 : if(optcell==3)then
259 11 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
260 11 : scale=(ucvol0/ucvol)**third
261 143 : rprimd(:,:)=scale*rprimd(:,:)
262 : end if
263 160 : call strainsym(nsym,rprimd0,rprimd,rprimd_symm,symrel)
264 640 : do jj=1,3
265 2080 : do ii=1,3
266 : ! write(ab_out,*) 'DIFF',ii,jj,abs(rprimd0(ii,jj)-rprimd_symm(ii,jj))
267 1440 : if (abs(rprimd0(ii,jj)-rprimd_symm(ii,jj))>1.E-14)&
268 1335 : & equal=.FALSE.
269 : end do
270 : end do
271 :
272 160 : if (equal)then
273 0 : acell(:)=acell0(:)
274 : rprimd(:,:)=rprimd0(:,:)
275 : else
276 : ! Use a representation based on normalised rprim vectors
277 160 : call mkradim(acell,rprim,rprimd_symm)
278 : end if
279 :
280 : end if
281 :
282 1009 : end subroutine xfpack_vin2x
283 : !!***
284 :
285 : !!****f* ABINIT/xfpack_x2vin
286 : !! NAME
287 : !! xfpack_x2vin
288 : !!
289 : !! FUNCTION
290 : !! Old option=1, transfer xred, acell, and rprim to vin
291 : !!
292 : !! INPUTS
293 : !! acell0(3)=reference length scales of primitive translations (bohr), needed for some values of optcell.
294 : !! natom=number of atoms in cell
295 : !! ndim=dimension of vin arrays
296 : !! nsym=order of group.
297 : !! rprimd0(3,3)=reference real space primitive translations,
298 : !! needed for some values of optcell.
299 : !! optcell=option for the optimisation of the unit cell. Described in abinit_help.
300 : !! Depending on its value, different part of acell and rprim
301 : !! are contained in vin.
302 : !! symrel(3,3,nsym)=symmetry operators in terms of action on primitive translations
303 : !! ucvol=unit cell volume (bohr^3), needed for some values of optcell.
304 : !! ucvol0=reference unit cell volume (bohr^3), needed for some values of optcell.
305 : !!
306 : !! OUTPUT
307 : !! (see side effects)
308 : !!
309 : !! SIDE EFFECTS
310 : !! Input/Output variables
311 : !! acell(3)=length scales of primitive translations (bohr)
312 : !! rprim(3,3)=dimensionless real space primitive translations
313 : !! vin(ndim)=vector that contains xred and some quantity derived
314 : !! from acell and rprim, depending on the value of optcell.
315 : !! xred(3,natom)=reduced dimensionless atomic coordinates
316 : !!
317 : !! SOURCE
318 :
319 1027 : subroutine xfpack_x2vin(acell,natom,ndim,nsym,optcell,&
320 1027 : & rprim,rprimd0,symrel,ucvol,ucvol0,vin,xred)
321 :
322 : !Arguments ------------------------------------
323 : !scalars
324 : integer,intent(in) :: natom,ndim,nsym,optcell
325 : real(dp),intent(in) :: ucvol0
326 : real(dp),intent(inout) :: ucvol !vz_i
327 : !arrays
328 : integer,intent(in) :: symrel(3,3,nsym)
329 : real(dp),intent(in) :: rprimd0(3,3)
330 : real(dp),intent(in) :: acell(3),rprim(3,3)
331 : real(dp),intent(in) :: xred(3,natom)
332 : real(dp),intent(out) :: vin(ndim)
333 :
334 : !Local variables-------------------------------
335 : !scalars
336 : integer :: ii,jj,kk
337 : real(dp) :: scale
338 : character(len=500) :: message
339 : !arrays
340 : real(dp) :: gmet(3,3),gprimd(3,3),gprimd0(3,3),rmet(3,3),rprimd(3,3)
341 : real(dp) :: rprimd_symm(3,3),scaling(3,3)
342 :
343 : ! *************************************************************************
344 :
345 : !!DEBUG
346 : !write(ab_out,*) ''
347 : !write(ab_out,*) 'xfpack_x2vin'
348 : !write(ab_out,*) 'natom=',natom
349 : !write(ab_out,*) 'ndim=',ndim
350 : !write(ab_out,*) 'nsym=',nsym
351 : !write(ab_out,*) 'optcell=',optcell
352 : !write(ab_out,*) 'ucvol=',ucvol
353 : !write(ab_out,*) 'xred='
354 : !do kk=1,natom
355 : !write(ab_out,*) xred(:,kk)
356 : !end do
357 : !write(ab_out,*) 'VECTOR INPUT (vin) xfpack_x2vin INPUT'
358 : !do ii=1,ndim,3
359 : !if (ii+2<=ndim)then
360 : !write(ab_out,*) ii,vin(ii:ii+2)
361 : !else
362 : !write(ab_out,*) ii,vin(ii:ndim)
363 : !end if
364 : !end do
365 : !!DEBUG
366 :
367 :
368 : !##########################################################
369 : !### 1. Test for compatible ndim
370 :
371 1027 : if(optcell==0 .and. ndim/=3*natom)then
372 : write(message,'(a,a,a,i4,a,i4,a)' )&
373 0 : & ' When optcell=0, ndim MUST be equal to 3*natom,',ch10,&
374 0 : & ' while ndim=',ndim,' and 3*natom=',3*natom,'.'
375 0 : ABI_BUG(message)
376 : end if
377 :
378 1027 : if( optcell==1 .and. ndim/=3*natom+1)then
379 : write(message,'(a,a,a,i4,a,i4,a)' )&
380 0 : & ' When optcell=1, ndim MUST be equal to 3*natom+1,',ch10,&
381 0 : & ' while ndim=',ndim,' and 3*natom+1=',3*natom+1,'.'
382 0 : ABI_BUG(message)
383 : end if
384 :
385 : if( (optcell==2 .or. optcell==3) &
386 1027 : & .and. ndim/=3*natom+6)then
387 : write(message,'(a,a,a,i4,a,i4,a)' )&
388 0 : & ' When optcell=2,3,4,5,6, ndim MUST be equal to 3*natom+6,',ch10,&
389 0 : & ' while ndim=',ndim,' and 3*natom+6=',3*natom+6,'.'
390 0 : ABI_BUG(message)
391 : end if
392 :
393 1027 : if( optcell>=4 .and. ndim/=3*natom+3)then
394 : write(message,'(a,a,a,i4,a,i4,a)' )&
395 0 : & ' When optcell=4,5,6,7,8 or 9, ndim MUST be equal to 3*natom+3,',ch10,&
396 0 : & ' while ndim=',ndim,' and 3*natom+3=',3*natom+3,'.'
397 0 : ABI_BUG(message)
398 : end if
399 :
400 : !##########################################################
401 : !### 2. option=1, transfer xred, acell, and rprim to vin
402 :
403 : !Get vin from xred, acell, and rprim
404 2054 : vin(1:3*natom)= reshape(xred(:,:), (/3*natom/) )
405 :
406 1027 : if(optcell/=0)then
407 318 : call mkrdim(acell,rprim,rprimd)
408 318 : call strainsym(nsym,rprimd0,rprimd,rprimd_symm,symrel)
409 318 : call metric(gmet,gprimd,-1,rmet,rprimd_symm,ucvol)
410 :
411 318 : if(optcell==1)then
412 :
413 : ! vin(3*natom+1)=ucvol**third
414 158 : vin(3*natom+1)=(ucvol/ucvol0)**third
415 :
416 160 : else if(optcell>=2)then
417 :
418 : ! Generates gprimd0
419 160 : call matr3inv(rprimd0,gprimd0)
420 160 : if (optcell==2 .or. optcell==3 .or. optcell>=7)then
421 584 : do ii=1,3
422 1898 : do jj=1,3
423 1314 : scaling(ii,jj)=0.0_dp
424 5694 : do kk=1,3
425 5256 : scaling(ii,jj)=scaling(ii,jj)+rprimd_symm(ii,kk)*gprimd0(jj,kk)
426 : end do
427 : end do
428 : end do
429 : ! Rescale if the volume must be preserved
430 146 : if(optcell==3)then
431 11 : scale=(ucvol0/ucvol)**third
432 143 : scaling(:,:)=scale*scaling(:,:)
433 : end if
434 146 : if(optcell==2 .or. optcell==3)then
435 129 : vin(3*natom+1)=scaling(1,1) ; vin(3*natom+4)=(scaling(2,3)+scaling(3,2))*0.5_dp
436 129 : vin(3*natom+2)=scaling(2,2) ; vin(3*natom+5)=(scaling(1,3)+scaling(3,1))*0.5_dp
437 129 : vin(3*natom+3)=scaling(3,3) ; vin(3*natom+6)=(scaling(1,2)+scaling(2,1))*0.5_dp
438 17 : else if(optcell>=7)then
439 17 : vin(3*natom+1)=scaling(1,1)
440 17 : vin(3*natom+2)=scaling(2,2)
441 17 : vin(3*natom+3)=scaling(3,3)
442 17 : if(optcell==7)vin(3*natom+1)=(scaling(2,3)+scaling(3,2))*0.5_dp
443 17 : if(optcell==8)vin(3*natom+2)=(scaling(1,3)+scaling(3,1))*0.5_dp
444 17 : if(optcell==9)vin(3*natom+3)=(scaling(1,2)+scaling(2,1))*0.5_dp
445 : end if
446 : end if
447 :
448 160 : if (optcell==4)then
449 9 : scaling(:,:) = 0.0_dp
450 36 : do ii=1,3
451 117 : do kk=1,3
452 108 : scaling(ii,1) = scaling(ii,1) + gprimd0(kk,ii)*rprimd_symm(kk,1)
453 : end do
454 : end do
455 9 : vin(3*natom+1) = scaling(1,1)
456 9 : vin(3*natom+2) = scaling(3,1)
457 9 : vin(3*natom+3) = scaling(2,1)
458 151 : else if (optcell==5)then
459 0 : scaling(:,:) = 0.0_dp
460 0 : do ii=1,3
461 0 : do kk=1,3
462 0 : scaling(ii,2) = scaling(ii,2) + gprimd0(kk,ii)*rprimd_symm(kk,2)
463 : end do
464 : end do
465 0 : vin(3*natom+1) = scaling(2,2)
466 0 : vin(3*natom+2) = scaling(3,2)
467 0 : vin(3*natom+3) = scaling(1,2)
468 151 : else if (optcell==6)then
469 5 : scaling(:,:) = 0.0_dp
470 20 : do ii=1,3
471 65 : do kk=1,3
472 60 : scaling(ii,3) = scaling(ii,3) + gprimd0(kk,ii)*rprimd_symm(kk,3)
473 : end do
474 : end do
475 5 : vin(3*natom+1) = scaling(3,3)
476 5 : vin(3*natom+2) = scaling(2,3)
477 5 : vin(3*natom+3) = scaling(1,3)
478 : end if
479 :
480 : end if
481 :
482 : end if
483 :
484 1027 : end subroutine xfpack_x2vin
485 : !!***
486 :
487 : !!****f* ABINIT/xfpack_f2vout
488 : !! NAME
489 : !! xfpack_f2vout
490 : !!
491 : !! FUNCTION
492 : !! Old option=3, transfer gred and strten to vout
493 : !!
494 : !! INPUTS
495 : !! natom=number of atoms in cell
496 : !! ndim=dimension of vout arrays
497 : !! optcell=option for the optimisation of the unit cell. Described in abinit_help.
498 : !! Depending on its value, different part of strten
499 : !! are contained in vout.
500 : !! strtarget(6)=target stresses ; they will be subtracted from strten when vout
501 : !! is computed.
502 : !! ucvol=unit cell volume (bohr^3), needed for some values of optcell.
503 : !!
504 : !! OUTPUT
505 : !! (see side effects)
506 : !!
507 : !! SIDE EFFECTS
508 : !! Input/Output variables
509 : !! gred(3,natom)=grads of Etot wrt reduced coordinates (hartree)
510 : !! strten(6)=components of the stress tensor (hartree/bohr^3)
511 : !! vout(ndim)=vector that contains gred and some quantity derived from
512 : !! strten, depending on the value of optcell, and taking care ot strtarget
513 : !!
514 : !! SOURCE
515 :
516 1009 : subroutine xfpack_f2vout(gred,natom,ndim,optcell,strtarget,strten,ucvol,vout)
517 :
518 : !Arguments ------------------------------------
519 : !scalars
520 : integer,intent(in) :: natom,ndim,optcell
521 : real(dp),intent(in) :: ucvol
522 : !arrays
523 : real(dp),intent(in) :: strtarget(6)
524 : real(dp),intent(in) :: gred(3,natom),strten(6)
525 : real(dp),intent(out) :: vout(ndim)
526 :
527 : !Local variables-------------------------------
528 : !scalars
529 : real(dp) :: strdiag
530 : character(len=500) :: message
531 : !arrays
532 : real(dp) :: dstr(6)
533 :
534 : ! *************************************************************************
535 :
536 : !!DEBUG
537 : !write(ab_out,*) ''
538 : !write(ab_out,*) 'xfpack_f2vout'
539 : !write(ab_out,*) 'natom=',natom
540 : !write(ab_out,*) 'ndim=',ndim
541 : !write(ab_out,*) 'optcell=',optcell
542 : !write(ab_out,*) 'ucvol=',ucvol
543 : !!DEBUG
544 :
545 :
546 : !##########################################################
547 : !### 1. Test for compatible ndim
548 :
549 1009 : if(optcell==0 .and. ndim/=3*natom)then
550 : write(message,'(a,a,a,i4,a,i4,a)' )&
551 0 : & ' When optcell=0, ndim MUST be equal to 3*natom,',ch10,&
552 0 : & ' while ndim=',ndim,' and 3*natom=',3*natom,'.'
553 0 : ABI_BUG(message)
554 : end if
555 :
556 1009 : if( optcell==1 .and. ndim/=3*natom+1)then
557 : write(message,'(a,a,a,i4,a,i4,a)' )&
558 0 : & ' When optcell=1, ndim MUST be equal to 3*natom+1,',ch10,&
559 0 : & ' while ndim=',ndim,' and 3*natom+1=',3*natom+1,'.'
560 0 : ABI_BUG(message)
561 : end if
562 :
563 : if( (optcell==2 .or. optcell==3) &
564 1009 : & .and. ndim/=3*natom+6)then
565 : write(message,'(a,a,a,i4,a,i4,a)' )&
566 0 : & ' When optcell=2 or 3, ndim MUST be equal to 3*natom+6,',ch10,&
567 0 : & ' while ndim=',ndim,' and 3*natom+6=',3*natom+6,'.'
568 0 : ABI_BUG(message)
569 : end if
570 :
571 1009 : if( optcell>=4 .and. ndim/=3*natom+3)then
572 : write(message,'(a,a,a,i4,a,i4,a)' )&
573 0 : & ' When optcell=4,5,6,7,8 or 9, ndim MUST be equal to 3*natom+3,',ch10,&
574 0 : & ' while ndim=',ndim,' and 3*natom+3=',3*natom+3,'.'
575 0 : ABI_BUG(message)
576 : end if
577 : !
578 : !Get vout from gred and strten
579 : !
580 2018 : vout(1:3*natom)= reshape(gred(:,:), (/3*natom/) )
581 7063 : dstr(:)=strten(:)-strtarget(:)
582 :
583 1009 : if(optcell==1)then
584 :
585 158 : vout(3*natom+1)=( dstr(1)+dstr(2)+dstr(3))*ucvol
586 :
587 851 : else if(optcell>=2)then
588 : ! Eventually take away the trace
589 160 : strdiag=0.0_dp
590 160 : if(optcell==3) strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp
591 : if(optcell==2 .or. optcell==3)then
592 516 : vout(3*natom+1:3*natom+3)=(dstr(1:3)-strdiag)*ucvol
593 : ! For non-diagonal derivatives, must take into account
594 : ! that eps(i,j) AND eps(j,i) are varied at the same time. Thus, derivative
595 : ! is twice larger
596 516 : vout(3*natom+4:3*natom+6)=dstr(4:6)*ucvol*2.0_dp
597 : else if(optcell==7 .or. optcell==8 .or. optcell==9)then
598 : ! Similar to case optcell==2 or optcell==3, but in 2 dimensions.
599 68 : vout(3*natom+1:3*natom+3)=dstr(1:3)*ucvol
600 17 : vout(3*natom+optcell-6) =dstr(optcell-3)*ucvol*2.0_dp
601 : else if (optcell==4)then
602 9 : vout(3*natom+1) = dstr(1)*ucvol
603 9 : vout(3*natom+2) = dstr(5)*ucvol
604 9 : vout(3*natom+3) = dstr(6)*ucvol
605 : else if (optcell==5)then
606 0 : vout(3*natom+1) = dstr(2)*ucvol
607 0 : vout(3*natom+2) = dstr(4)*ucvol
608 0 : vout(3*natom+3) = dstr(6)*ucvol
609 : else if (optcell==6)then
610 5 : vout(3*natom+1) = dstr(3)*ucvol
611 5 : vout(3*natom+2) = dstr(4)*ucvol
612 5 : vout(3*natom+3) = dstr(5)*ucvol
613 : end if
614 :
615 : end if
616 :
617 1009 : end subroutine xfpack_f2vout
618 : !!***
619 :
620 :
621 : !!****f* ABINIT/xfh_recover_new
622 : !! NAME
623 : !! xfh_recover_new
624 : !!
625 : !! FUNCTION
626 : !! Update the contents of the history xfhist taking values
627 : !! from xred, acell, rprim, gred_corrected and strten
628 : !!
629 : !! INPUTS
630 : !!
631 : !! OUTPUT
632 : !!
633 : !! SOURCE
634 :
635 :
636 7 : subroutine xfh_recover_new(ab_xfh,ab_mover,acell,cycl_main,gred,&
637 7 : & hessin,ndim,rprim,rprimd0,strten,ucvol,ucvol0,vin,vin_prev,vout,&
638 14 : & vout_prev,xred)
639 :
640 : !Arguments ------------------------------------
641 : !scalars
642 :
643 : integer,intent(in) :: ndim
644 : integer,intent(out) :: cycl_main
645 : real(dp),intent(inout) :: ucvol,ucvol0
646 : type(ab_xfh_type),intent(inout) :: ab_xfh
647 : type(abimover),intent(in) :: ab_mover
648 :
649 :
650 : !arrays
651 : real(dp),intent(inout) :: acell(3)
652 : real(dp),intent(inout) :: hessin(:,:)
653 : real(dp),intent(inout) :: xred(3,ab_mover%natom)
654 : real(dp),intent(inout) :: rprim(3,3)
655 : real(dp),intent(inout) :: rprimd0(3,3)
656 : real(dp),intent(inout) :: gred(3,ab_mover%natom)
657 : real(dp),intent(inout) :: strten(6)
658 : real(dp),intent(inout) :: vin(:)
659 : real(dp),intent(inout) :: vin_prev(:)
660 : real(dp),intent(inout) :: vout(:)
661 : real(dp),intent(inout) :: vout_prev(:)
662 :
663 : !Local variables-------------------------------
664 : !scalars
665 : integer :: ixfh ! kk,jj
666 :
667 : !*********************************************************************
668 :
669 7 : if(ab_xfh%nxfh/=0)then
670 : ! Loop over previous time steps
671 0 : do ixfh=1,ab_xfh%nxfh
672 :
673 : ! For that time step, get new (x,f) from xfhist
674 0 : xred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom ,1,ixfh)
675 0 : rprim(1:3,1:3)=ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh)
676 0 : acell(:) =ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh)
677 0 : gred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh)
678 : ! This use of results_gs is unusual
679 0 : strten(1:3) =ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh)
680 0 : strten(4:6) =ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh)
681 :
682 : ! !DEBUG
683 : ! write (ab_out,*) '---READED FROM XFHIST---'
684 :
685 : ! write (ab_out,*) 'XRED'
686 : ! do kk=1,ab_mover%natom
687 : ! write (ab_out,*) xred(:,kk)
688 : ! end do
689 : ! write (ab_out,*) 'FRED'
690 : ! do kk=1,ab_mover%natom
691 : ! write (ab_out,*) gred(:,kk)
692 : ! end do
693 : ! write(ab_out,*) 'RPRIM'
694 : ! do kk=1,3
695 : ! write(ab_out,*) rprim(:,kk)
696 : ! end do
697 : ! write(ab_out,*) 'ACELL'
698 : ! write(ab_out,*) acell(:)
699 : ! !DEBUG
700 :
701 : ! Transfer it in vin, vout
702 : call xfpack_x2vin(acell,ab_mover%natom,&
703 : & ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,&
704 0 : & ab_mover%symrel,ucvol,ucvol0,vin,xred)
705 : call xfpack_f2vout(gred,ab_mover%natom,&
706 : & ndim,ab_mover%optcell,ab_mover%strtarget,strten,&
707 0 : & ucvol,vout)
708 : ! Get old time step, if any, and update inverse hessian
709 0 : if(ixfh/=1)then
710 0 : xred(:,:) =ab_xfh%xfhist(:,1:ab_mover%natom,1,ixfh-1)
711 : rprim(1:3,1:3)=&
712 0 : & ab_xfh%xfhist(:,ab_mover%natom+2:ab_mover%natom+4,1,ixfh-1)
713 0 : acell(:)=ab_xfh%xfhist(:,ab_mover%natom+1,1,ixfh-1)
714 0 : gred(:,:)=ab_xfh%xfhist(:,1:ab_mover%natom,2,ixfh-1)
715 : ! This use of results_gs is unusual
716 0 : strten(1:3)=ab_xfh%xfhist(:,ab_mover%natom+2,2,ixfh-1)
717 0 : strten(4:6)=ab_xfh%xfhist(:,ab_mover%natom+3,2,ixfh-1)
718 : ! Tranfer it in vin_prev, vout_prev
719 : call xfpack_x2vin(acell,ab_mover%natom,&
720 : & ndim,ab_mover%nsym,ab_mover%optcell,rprim,rprimd0,&
721 0 : & ab_mover%symrel,ucvol,ucvol0,vin_prev,xred)
722 : call xfpack_f2vout(gred,ab_mover%natom,&
723 : & ndim,ab_mover%optcell,ab_mover%strtarget,strten,&
724 0 : & ucvol,vout_prev)
725 :
726 : ! write(ab_out,*) 'Hessian matrix before update',ndim,'x',ndim
727 : ! write(ab_out,*) 'ixfh=',ixfh
728 : ! do kk=1,ndim
729 : ! do jj=1,ndim,3
730 : ! if (jj+2<=ndim)then
731 : ! write(ab_out,*) jj,hessin(jj:jj+2,kk)
732 : ! else
733 : ! write(ab_out,*) jj,hessin(jj:ndim,kk)
734 : ! end if
735 : ! end do
736 : ! end do
737 :
738 : call hessupdt(hessin,ab_mover%iatfix,ab_mover%natom,ndim,&
739 0 : & vin,vin_prev,vout,vout_prev)
740 :
741 : ! !DEBUG
742 : ! write(ab_out,*) 'Hessian matrix after update',ndim,'x',ndim
743 : ! do kk=1,ndim
744 : ! do jj=1,ndim,3
745 : ! if (jj+2<=ndim)then
746 : ! write(ab_out,*) jj,hessin(jj:jj+2,kk)
747 : ! else
748 : ! write(ab_out,*) jj,hessin(jj:ndim,kk)
749 : ! end if
750 : ! end do
751 : ! end do
752 : ! !DEBUG
753 :
754 : end if !if(ab_xfh%nxfh/=0)
755 : end do ! End loop over previous time steps
756 :
757 : ! The hessian has been generated,
758 : ! as well as the latest vin and vout
759 : ! so will cycle the main loop
760 0 : cycl_main=1
761 : end if
762 :
763 7 : end subroutine xfh_recover_new
764 : !!***
765 :
766 : !!****f* ABINIT/xfh_update
767 : !! NAME
768 : !! xfh_update
769 : !!
770 : !! FUNCTION
771 : !! Update the contents of the history xfhist taking values
772 : !! from xred, acell, rprim, gred_corrected and strten
773 : !!
774 : !! INPUTS
775 : !!
776 : !! OUTPUT
777 : !!
778 : !! SOURCE
779 :
780 1158 : subroutine xfh_update(ab_xfh,acell,gred_corrected,natom,rprim,strten,xred)
781 :
782 : !Arguments ------------------------------------
783 : !scalars
784 : type(ab_xfh_type),intent(inout) :: ab_xfh
785 : integer,intent(in) :: natom
786 :
787 : !arrays
788 : real(dp),intent(in) :: acell(3)
789 : real(dp),intent(in) :: xred(3,natom)
790 : real(dp),intent(in) :: rprim(3,3)
791 : real(dp),intent(in) :: gred_corrected(3,natom)
792 : real(dp),intent(in) :: strten(6)
793 :
794 : !Local variables-------------------------------
795 : !scalars
796 : !integer :: kk
797 :
798 : !*********************************************************************
799 :
800 : !DEBUG
801 : !write (ab_out,*) '---WROTE TO XFHIST---'
802 :
803 : !write (ab_out,*) 'XRED'
804 : !do kk=1,natom
805 : !write (ab_out,*) xred(:,kk)
806 : !end do
807 : !write (ab_out,*) 'FRED'
808 : !do kk=1,natom
809 : !write (ab_out,*) gred_corrected(:,kk)
810 : !end do
811 : !write(ab_out,*) 'RPRIM'
812 : !do kk=1,3
813 : !write(ab_out,*) rprim(:,kk)
814 : !end do
815 : !write(ab_out,*) 'ACELL'
816 : !write(ab_out,*) acell(:)
817 : !DEBUG
818 :
819 1158 : ab_xfh%nxfh=ab_xfh%nxfh+1
820 :
821 11986 : ab_xfh%xfhist(:,1:natom,1,ab_xfh%nxfh)=xred(:,:)
822 4632 : ab_xfh%xfhist(:,natom+1,1,ab_xfh%nxfh)=acell(:)
823 15054 : ab_xfh%xfhist(:,natom+2:natom+4,1,ab_xfh%nxfh)=rprim(:,:)
824 11986 : ab_xfh%xfhist(:,1:natom,2,ab_xfh%nxfh)=gred_corrected(:,:)
825 4632 : ab_xfh%xfhist(:,natom+2,2,ab_xfh%nxfh)=strten(1:3)
826 4632 : ab_xfh%xfhist(:,natom+3,2,ab_xfh%nxfh)=strten(4:6)
827 :
828 1158 : end subroutine xfh_update
829 : !!***
830 :
831 : end module m_xfpack
832 : !!***
|
