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