Branch data Line data Source code
1 : : !!****m* ABINIT/m_relaxpol
2 : : !! NAME
3 : : !! m_relaxpol
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1999-2026 ABINIT group (MVeithen)
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_relaxpol
22 : :
23 : : use defs_basis
24 : : use m_abicore
25 : : use m_errors
26 : :
27 : : use m_fstrings, only : sjoin, itoa
28 : : use m_matrix, only : matr3inv
29 : : use m_berrytk, only : polcart
30 : : use m_hide_lapack, only : dzgedi, dzgefa
31 : : use m_geometry, only : xcart2xred
32 : : use m_dynmat, only : symdyma
33 : : use m_crystal, only : crystal_t
34 : :
35 : : implicit none
36 : :
37 : : private
38 : : !!***
39 : :
40 : : public :: relaxpol
41 : : !!***
42 : :
43 : : contains
44 : : !!***
45 : :
46 : : !!****f* ABINIT/relaxpol
47 : : !! NAME
48 : : !! relaxpol
49 : : !!
50 : : !! FUNCTION
51 : : !! 1) Compute polarization in cartesian coordinates
52 : : !! 2) Structural relaxation at fixed polarization: this routine
53 : : !! solves the linear system of equations Eq.(13)
54 : : !! of Na Sai et al., PRB 66, 104108 (2002) [[cite:Sai2002]].
55 : : !!
56 : : !! INPUTS
57 : : !! blkflg(msize) = flag for every matrix element (0=> the element
58 : : !! is not in the data block), (1=> the element is in the data blok)
59 : : !! blkval(2,msize) = matrix that contains the second-order energy derivatives
60 : : !! etotal = Kohn-Sham energy at zero electric field
61 : : !! gred(3,natom) = -1 times the forces in reduced coordinates
62 : : !! iatfix(natom) = indices of the atoms that are held fixed in the relaxation
63 : : !! iout = unit number for output
64 : : !! istrfix(6) = indices of the elements of the strain tensor that
65 : : !! are held fixed in the relaxation
66 : : !! 1 = xx
67 : : !! 2 = yy
68 : : !! 3 = zz
69 : : !! 4 = yz & zy
70 : : !! 5 = xz & zx
71 : : !! 6 = xy & yx
72 : : !! mpert = maximum number of ipert
73 : : !! msize = dimension of blkflg and blkval
74 : : !! natfix = number of atoms that are held fixed in the relaxation
75 : : !! natom = number of atoms in the unit cell
76 : : !! nstrfix = number of elements of the strain tensor that are held fixed in the relaxation
77 : : !! pel(3) = electronic polarization not taking into account the factor 1/ucvol
78 : : !! red_ptot(3) = total polarization reduced units !!REC
79 : : !! relaxat = 1: relax atomic positions
80 : : !! = 0: do not relax atomic positions
81 : : !! relaxstr = 1: relax cell parameters
82 : : !! = 0: do not relax cell parameters
83 : : !! strten(6) = stress tensor in cartesian coordinates
84 : : !! targetpol(3) = target value of the polarization
85 : : !!
86 : : !! OUTPUT
87 : : !!
88 : : !! NOTES
89 : : !! - The elements of the dynamical matrix stored in blkval
90 : : !! are symmetrized before computing the new atomic positions and cell parameters.
91 : : !! - In case relaxat = 0 and relaxstr = 0, the routine only
92 : : !! computes the polarization in cartesian coordinates.
93 : : !!
94 : : !! SOURCE
95 : :
96 : 3 : subroutine relaxpol(Crystal,blkflg,blkval,etotal,gred,iatfix,iout,istrfix,&
97 : : & mpert,msize,natfix,natom,nstrfix,pel,red_ptot,relaxat,relaxstr,&
98 : : & strten,targetpol)
99 : :
100 : : !Arguments -------------------------------
101 : : !scalars
102 : : integer,intent(in) :: iout,mpert,msize,natfix,natom,nstrfix
103 : : integer,intent(in) :: relaxat,relaxstr
104 : : real(dp),intent(in) :: etotal
105 : : type(crystal_t),intent(in) :: Crystal
106 : : !arrays
107 : : integer,intent(in) :: blkflg(msize),iatfix(natom)
108 : : integer,intent(in) :: istrfix(6)
109 : : real(dp),intent(in) :: gred(3,natom),pel(3),strten(6)
110 : : real(dp),intent(in) :: red_ptot(3)
111 : : real(dp),intent(inout) :: blkval(2,msize),targetpol(3)
112 : :
113 : : !Local variables -------------------------
114 : : !scalars
115 : : integer :: flag,iatom,idir,ii,index,index1,index_tild,info,ipert,istrain
116 : : integer :: itypat,jdir,job,jpert,polunit,posi,posj,sizef
117 : : integer :: usepaw
118 : : logical :: iwrite
119 : : real(dp) :: e1,fmax,poltmp,sigmax,tol,value,ucvol
120 : : character(len=500) :: message
121 : : !arrays
122 : : integer :: irelaxstrain(6)
123 : 3 : integer,allocatable :: ipvt(:),irelaxat(:),rfpert(:,:)
124 : 6 : real(dp) :: acell_new(3),delta_eta(6),delta_xcart(3,natom),det(2,2),diffpol(3),rprimd(3,3)
125 : : real(dp) :: diffsig(6),favg(3),gprimd(3,3),lambda(3),pel_cart(3),pelev(3)
126 : : real(dp) :: pion(3),pion_cart(3),ptot_cart(3),qphon(3),rprim(3,3)
127 : : real(dp) :: rprimd_new(3,3),sigelfd(6),strainmat(3,3)
128 : 6 : real(dp) :: xcart_new(3,natom),xred_new(3,natom)
129 : 3 : real(dp),allocatable :: cfac(:,:),delta(:),dymati(:),fcart(:,:),fcmat(:,:,:)
130 : 3 : real(dp),allocatable :: fdiff(:,:),felfd(:,:),ifcmat(:,:,:),vec(:),zgwork(:,:)
131 : :
132 : : ! *********************************************************************
133 : :
134 : 3 : usepaw = 0
135 [ + + + + ]: 39 : rprimd = Crystal%rprimd
136 : 3 : ucvol = Crystal%ucvol
137 : 3 : iwrite = iout > 0
138 : :
139 : : !Check if some degrees of freedom remain fixed during the optimization
140 : :
141 [ - + + - : 9 : ABI_MALLOC(irelaxat,(natom))
+ - + - -
+ - + ]
142 [ + + + + ]: 30 : irelaxat(:) = 1 ; irelaxstrain(:) = 1
143 [ - + ]: 3 : if (natfix > 0) then
144 [ # # ]: 0 : do ii = 1, natfix
145 : 0 : iatom = iatfix(ii)
146 [ # # # # ]: 0 : if ((iatom > natom).or.(iatom < 0)) then
147 : : write(message, '(a,i0,a,i0,a,a,a,a,a)')&
148 : 0 : & 'The value of iatfix(',ii,') is ',iatom,', which is not allowed.',ch10,&
149 : 0 : & 'iatfix must be larger than 0 and smaller than natom.',ch10,&
150 : 0 : & 'Action: correct iatfix in your input file.'
151 : 0 : ABI_ERROR(message)
152 : : end if
153 : 0 : irelaxat(iatom) = 0
154 : : end do
155 : : end if
156 : :
157 [ + + ]: 3 : if (nstrfix > 0) then
158 [ + + ]: 3 : do ii = 1, nstrfix
159 : 2 : istrain = istrfix(ii)
160 [ - + ]: 2 : if ((istrain > 6).or.(istrain < 0)) then
161 : : write(message, '(a,i0,a,i0,a,a,a,a,a)')&
162 : 0 : & 'istrfix(',ii,') is',istrain,', which is not allowed.',ch10,&
163 : 0 : & 'istrfix must be larger than 0 and smaller than 6.',ch10,&
164 : 0 : & 'Action : correct istrfix in your input file.'
165 : 0 : ABI_ERROR(message)
166 : : end if
167 : 3 : irelaxstrain(istrain) = 0
168 : : end do
169 : : end if
170 : :
171 : :
172 [ - + + - : 9 : ABI_MALLOC(rfpert,(mpert,3))
+ - + - -
+ - + ]
173 [ - + + - : 12 : ABI_MALLOC(cfac,(mpert,mpert))
- + + - +
- + - - +
- + ]
174 : 3 : call matr3inv(rprimd,gprimd)
175 : :
176 : : !Compute the size of the matrix that contains the second-order derivatives
177 : :
178 : 3 : sizef = 3
179 [ + + + + ]: 165 : rfpert(:,:) = 0
180 [ + + ]: 12 : rfpert(natom+2,1:3) = 1
181 [ + + ]: 3 : if (relaxat == 1) then
182 [ + + ]: 9 : do iatom = 1, natom
183 [ + - ]: 9 : if (irelaxat(iatom) == 1) then
184 : 7 : sizef = sizef + 3
185 [ + + ]: 28 : rfpert(iatom,1:3) = 1
186 : : end if
187 : : end do
188 : : end if
189 : 3 : ii = natom + 2
190 [ + + ]: 3 : if (relaxstr == 1) then
191 : 2 : istrain = 0
192 [ + + ]: 6 : do ipert = (natom+3), (natom+4)
193 [ + + ]: 18 : do idir = 1, 3
194 : 12 : istrain = istrain + 1
195 [ + + ]: 16 : if (irelaxstrain(istrain) == 1) then
196 : 10 : sizef = sizef + 1
197 : 10 : rfpert(ipert,idir) = 1
198 : : end if
199 : : end do
200 : : end do
201 : : end if
202 : :
203 [ - + + - : 12 : ABI_MALLOC(fcmat,(2,sizef,sizef))
- + + - +
- + - - +
- + ]
204 [ - + + - : 9 : ABI_MALLOC(ifcmat,(2,sizef,sizef))
- + + - -
+ - + ]
205 [ - + + - : 9 : ABI_MALLOC(vec,(sizef))
+ - + - -
+ - + ]
206 [ - + + - : 6 : ABI_MALLOC(delta,(sizef))
- + - + ]
207 [ - + + - : 9 : ABI_MALLOC(ipvt,(sizef))
+ - + - -
+ - + ]
208 [ - + + - : 9 : ABI_MALLOC(zgwork,(2,sizef))
+ - + - -
+ - + ]
209 [ - + + - : 9 : ABI_MALLOC(fcart,(3,natom))
+ - + - -
+ - + ]
210 [ - + + - : 6 : ABI_MALLOC(felfd,(3,natom))
- + - + ]
211 [ - + + - : 6 : ABI_MALLOC(fdiff,(3,natom))
- + - + ]
212 : :
213 : : !Build the vector that stores the forces, sigma and the polarization
214 : :
215 [ + + ]: 43 : vec(:) = zero
216 : 3 : posi = 0
217 : :
218 [ + + ]: 3 : if (relaxat == 1) then
219 : :
220 : : ! Note conversion to cartesian coordinates (bohr) AND
221 : : ! negation to make a force out of a gradient
222 : : ! Also subtract off average force from each force
223 : : ! component to avoid spurious drifting of atoms across cell.
224 : 2 : favg(:) = zero
225 [ + + ]: 9 : do iatom = 1, natom
226 [ + + ]: 30 : do idir = 1, 3
227 : : fcart(idir,iatom) = -(gprimd(idir,1)*gred(1,iatom) + &
228 : : & gprimd(idir,2)*gred(2,iatom) + &
229 : 21 : & gprimd(idir,3)*gred(3,iatom))
230 : 28 : favg(idir) = favg(idir) + fcart(idir,iatom)
231 : : end do
232 : : end do
233 [ + + ]: 8 : favg(:) = favg(:)/dble(natom)
234 [ + + ]: 9 : do iatom = 1, natom
235 [ + + ]: 30 : fcart(:,iatom) = fcart(:,iatom) - favg(:)
236 : : end do
237 : :
238 [ + + ]: 9 : do iatom = 1, natom
239 [ + - ]: 9 : if (irelaxat(iatom) == 1) then
240 [ + + ]: 28 : do idir = 1, 3
241 : 21 : posi = posi + 1
242 : 28 : vec(posi) = fcart(idir,iatom)
243 : : end do
244 : : end if
245 : : end do
246 : :
247 : : end if ! relaxat == 1
248 : :
249 : : !DEBUG
250 : : !write(std_out,*)'Forces in cartesian coords'
251 : : !do iatom = 1, natom
252 : : !write(std_out,'(3(2x,e16.9))')(fcart(idir,iatom),idir = 1, 3)
253 : : !end do
254 : : !stop
255 : : !ENDDEBUG
256 : :
257 : : !Transform target polarization to atomic units
258 [ + + ]: 12 : targetpol(:) = targetpol(:)*((Bohr_Ang*1.0d-10)**2)/e_Cb
259 : :
260 : : !Compute ionic polarization
261 : 3 : pion(:) = zero
262 [ + + ]: 12 : do iatom = 1, natom
263 : 9 : itypat = Crystal%typat(iatom)
264 [ + + ]: 39 : do idir = 1, 3
265 : 27 : poltmp = Crystal%zion(itypat) * Crystal%xred(idir,iatom)
266 : 27 : poltmp = poltmp - two*nint(poltmp/two) ! fold into [-1,1]
267 : 36 : pion(idir) = pion(idir) + poltmp
268 : : end do
269 : : end do
270 [ + + ]: 12 : do idir = 1, 3
271 : 12 : pion(idir) = pion(idir) - two*nint(pion(idir)/two) ! fold into [-1,1]
272 : : end do
273 : :
274 : : !Transform the polarization to cartesian coordinates
275 : 3 : polunit = 3
276 : 3 : pelev=zero ! This is a PAW-related quantity, which we ignore here.
277 : : call polcart(red_ptot,pel,pel_cart,pelev,pion,pion_cart,polunit,&
278 : 3 : & ptot_cart,rprimd,ucvol,iout,usepaw)
279 : :
280 [ + + ]: 12 : do idir = 1, 3
281 : 9 : posi = posi + 1
282 : 12 : vec(posi) = ptot_cart(idir) - targetpol(idir)
283 : : end do
284 : :
285 : :
286 [ + + ]: 3 : if (relaxstr == 1) then
287 [ + + ]: 14 : do istrain = 1, 6
288 [ + + ]: 14 : if (irelaxstrain(istrain) == 1) then
289 : 10 : posi = posi + 1
290 : 10 : vec(posi) = -1._dp*strten(istrain)*ucvol
291 : : end if
292 : : end do
293 : : end if
294 : :
295 : :
296 : : !Symmetrize the dynamical matrix
297 : :
298 [ + - - + : 9 : ABI_MALLOC(dymati,(2*3*natom*3*natom))
+ - + - +
- - + -
+ ]
299 : : !by the symdyma routine
300 [ + + ]: 12 : do ipert = 1, natom
301 [ + + ]: 39 : do idir = 1, 3
302 [ + + ]: 135 : do jpert = 1, natom
303 [ + + ]: 423 : do jdir = 1, 3
304 : 297 : index = jdir +3*((jpert - 1) + mpert*((idir - 1) + 3*(ipert - 1)))
305 : 297 : index1 = jdir +3*((jpert - 1) + natom*((idir - 1) + 3*(ipert - 1)))
306 : 297 : dymati(2*index1 - 1) = blkval(1,index)
307 : 396 : dymati(2*index1 ) = blkval(2,index)
308 : : end do
309 : : end do
310 : : end do
311 : : end do
312 : :
313 : 3 : qphon(:) = zero
314 : 3 : call symdyma(dymati,Crystal%indsym,natom,Crystal%nsym,qphon,rprimd,Crystal%symrel,Crystal%symafm)
315 : :
316 [ + + ]: 12 : do ipert = 1, natom
317 [ + + ]: 39 : do idir = 1, 3
318 [ + + ]: 135 : do jpert = 1, natom
319 [ + + ]: 423 : do jdir = 1, 3
320 : 297 : index = jdir +3*((jpert - 1) + mpert*((idir - 1) + 3*(ipert - 1)))
321 : 297 : index1 = jdir +3*((jpert - 1) + natom*((idir - 1) + 3*(ipert - 1)))
322 : 297 : blkval(1,index) = dymati(2*index1 - 1)
323 : 396 : blkval(2,index) = dymati(2*index1 )
324 : : end do
325 : : end do
326 : : end do
327 : : end do
328 : :
329 [ - + ]: 3 : ABI_FREE(dymati)
330 : :
331 : : !Define conversion factors for blkval
332 [ + + + + ]: 945 : cfac(:,:) = 1._dp
333 [ + + ]: 12 : cfac(1:natom,natom+2) = -1._dp/ucvol
334 [ + + ]: 12 : cfac(natom+2,1:natom) = -1._dp/ucvol
335 [ + + ]: 9 : cfac(natom+3:natom+4,natom+2) = -1._dp
336 [ + + ]: 9 : cfac(natom+2,natom+3:natom+4) = -1._dp
337 : :
338 : :
339 : : !Build the matrix that contains the second-order derivatives
340 : : !ipert = natom + 1 corresponds to the ddk perturbation, that
341 : : !is not needed; so skip it
342 : :
343 [ + + + + : 2197 : fcmat(:,:,:) = zero
+ + ]
344 : :
345 : 3 : posi = 0
346 : 3 : flag = 0
347 : : ! When fcmat has been build, flag = 0 if all elements were available.
348 : : ! Otherwise, it will be 1. In case one element is missing, check if
349 : : ! it can be obtained by changing the order of the perturbations
350 : :
351 [ + + ]: 54 : do ipert = 1, mpert
352 [ + + ]: 207 : do idir = 1, 3
353 [ + + ]: 204 : if (rfpert(ipert,idir) == 1) then
354 : 40 : posi = posi + 1
355 : 40 : posj = 0
356 : :
357 [ + + ]: 772 : do jpert = 1, mpert
358 [ + + ]: 2968 : do jdir = 1, 3
359 [ + + ]: 2928 : if (rfpert(jpert,jdir) == 1) then
360 : 718 : index = jdir +3*((jpert - 1) + mpert*((idir - 1) + 3*(ipert - 1)))
361 : 718 : index_tild = idir +3*((ipert - 1) + mpert*((jdir - 1) + 3*(jpert - 1)))
362 : 718 : posj = posj + 1
363 [ + + + + ]: 718 : if ((ipert /= natom + 2).or.(jpert /= natom + 2)) then
364 [ + - ]: 691 : if (blkflg(index) == 1) then
365 [ + + ]: 2073 : fcmat(:,posi,posj) = blkval(:,index)*cfac(ipert,jpert)
366 [ # # ]: 0 : else if (blkflg(index_tild) == 1) then
367 [ # # ]: 0 : fcmat(:,posi,posj) = blkval(:,index_tild)*cfac(ipert,jpert)
368 [ # # ]: 0 : blkval(:,index) = blkval(:,index_tild)
369 : : else
370 : 0 : flag = 1
371 : 0 : write(std_out,'(a,4(2x,i3))')'relaxpol: could not find element:',idir,ipert,jdir,jpert
372 : : end if
373 : : end if
374 : : ! DEBUG
375 : : ! write(100,'(4(2x,i3),5x,f16.9)')idir,ipert,jdir,jpert,fcmat(1,posi,posj)
376 : : ! ENDDEBUG
377 : : end if
378 : : end do
379 : : end do
380 : :
381 : : end if
382 : : end do
383 : : end do
384 : :
385 [ - + ]: 3 : if (flag == 1) then
386 : : write(message, '(a,a,a,i0,a,i0,a,a,a,a)' )&
387 : 0 : & 'Some of the second order derivatives required to deal with the case',ch10,&
388 : 0 : & 'relaxat = ',relaxat,', relaxstr = ', relaxstr, ch10,&
389 : 0 : & 'are missing in the DDB.',ch10,&
390 : 0 : & 'Action: correct your DDB or change your input file.'
391 : 0 : ABI_ERROR(message)
392 : : end if
393 : :
394 : :
395 : : !Compute the inverse of the force constant matrix
396 : :
397 [ + + - + ]: 3 : if ((relaxat /= 0).or.(relaxstr /= 0)) then
398 : :
399 : 2 : job = 1 ! compute inverse only
400 [ + + + + : 2166 : ifcmat(:,:,:) = fcmat(:,:,:)
+ + ]
401 : :
402 : 2 : call dzgefa(ifcmat,sizef,sizef,ipvt,info)
403 [ - + ]: 2 : ABI_CHECK(info == 0, sjoin("dzgefa returned:", itoa(info)))
404 : 2 : call dzgedi(ifcmat,sizef,sizef,ipvt,det,zgwork,job)
405 : :
406 : : ! DEBUG
407 : : ! write(100,*)'relaxat = ',relaxat
408 : : ! write(100,*)'relaxstr = ',relaxstr
409 : : ! write(100,*)'irelaxat = '
410 : : ! write(100,*)irelaxat(:)
411 : : ! write(100,*)'irelaxstrain = '
412 : : ! write(100,*)irelaxstrain(:)
413 : : ! write(100,*)'sizef = ',sizef
414 : : ! write(100,*)'targetpol ='
415 : : ! write(100,*)targetpol(:)
416 : : ! do ipert = 1, sizef
417 : : ! do jpert = 1, sizef
418 : : ! write(100,'(2(2x,i3),2x,e16.9)')ipert,jpert,fcmat(1,ipert,jpert)
419 : : ! end do
420 : : ! end do
421 : : ! stop
422 : : ! ENDDEBUG
423 : :
424 : : ! Compute \delta R, \delta \eta and \lambda
425 [ + + ]: 39 : delta(:) = zero
426 [ + + ]: 39 : do ipert = 1, sizef
427 [ + + ]: 748 : do jpert = 1, sizef
428 : 746 : delta(ipert) = delta(ipert) + ifcmat(1,ipert,jpert)*vec(jpert)
429 : : end do
430 : : end do
431 : :
432 : :
433 : : ! Update atomic positions
434 : 2 : posi = 0
435 [ + - ]: 2 : if (relaxat == 1) then
436 : :
437 [ + + + + ]: 30 : delta_xcart(:,:) = zero
438 [ + + + + ]: 30 : xcart_new(:,:) = zero
439 [ + + ]: 9 : do iatom = 1, natom
440 [ + - ]: 9 : if (irelaxat(iatom) == 1) then
441 [ + + ]: 28 : do idir = 1, 3
442 : 21 : posi = posi + 1
443 : 28 : delta_xcart(idir,iatom) = delta(posi)
444 : : end do
445 : : end if
446 : : end do
447 : :
448 : : ! Drop unsignificant digits in order to eleminate numerical noise
449 : 2 : tol = 10000000._dp
450 [ + + ]: 9 : do iatom = 1, natom
451 [ + + ]: 30 : do idir = 1, 3
452 : 21 : value = delta_xcart(idir,iatom)
453 : 21 : ii = log10(abs(value))
454 [ + - ]: 21 : if (ii <= 0) then
455 : 21 : ii = abs(ii) + 1
456 : 21 : value = one*int(tol*value*10.0_dp**ii)/(tol*10.0_dp**ii) !vz_d
457 : : else
458 : 0 : value = one*int(tol*value/(10.0_dp**ii))*(10.0_dp**ii)/tol !vz_d
459 : : end if
460 : 28 : delta_xcart(idir,iatom) = value
461 : : end do
462 : : end do
463 : :
464 [ + + + + ]: 30 : xcart_new(:,:) = Crystal%xcart(:,:) + delta_xcart(:,:)
465 : 2 : call xcart2xred(natom,rprimd,xcart_new,xred_new)
466 : : end if ! relaxat == 1
467 : :
468 : : ! Compute lambda and the value of the energy functional F - \lambda \cdot P$
469 : :
470 : 2 : e1 = etotal
471 [ + + ]: 8 : do idir = 1, 3
472 : 6 : posi = posi + 1
473 : 6 : lambda(idir) = delta(posi)
474 : 8 : e1 = e1 - lambda(idir)*ptot_cart(idir)
475 : : end do
476 : :
477 : : ! Update cell parameters
478 [ + - ]: 2 : if (relaxstr == 1) then
479 : 2 : delta_eta(:) = zero
480 [ + + ]: 14 : do istrain = 1, 6
481 [ + + ]: 14 : if (irelaxstrain(istrain) == 1) then
482 : 10 : posi = posi + 1
483 : 10 : delta_eta(istrain) = delta(posi)
484 : : end if
485 : : end do
486 : :
487 [ + + ]: 8 : do istrain = 1, 3
488 : 8 : strainmat(istrain,istrain) = delta_eta(istrain)
489 : : end do
490 : 2 : strainmat(2,3) = delta_eta(4)/2._dp ; strainmat(3,2) = delta_eta(4)/2._dp
491 : 2 : strainmat(1,3) = delta_eta(5)/2._dp ; strainmat(3,1) = delta_eta(5)/2._dp
492 : 2 : strainmat(2,1) = delta_eta(6)/2._dp ; strainmat(1,2) = delta_eta(6)/2._dp
493 : :
494 : 2 : rprimd_new(:,:) = 0._dp
495 [ + + ]: 8 : do idir = 1, 3
496 [ + + ]: 26 : do jdir = 1, 3
497 [ + + ]: 78 : do ii = 1, 3
498 : : rprimd_new(jdir,idir) = rprimd_new(jdir,idir) + &
499 : 72 : & rprimd(ii,idir)*strainmat(ii,jdir)
500 : : end do
501 : : end do
502 : : end do
503 [ + + + + ]: 26 : rprimd_new(:,:) = rprimd_new(:,:) + rprimd(:,:)
504 : :
505 : 2 : acell_new(:) = zero
506 [ + + ]: 8 : do idir = 1, 3
507 [ + + ]: 24 : do jdir = 1, 3
508 : : acell_new(idir) = acell_new(idir) + &
509 : 24 : & rprimd_new(jdir,idir)*rprimd_new(jdir,idir)
510 : : end do
511 : 6 : acell_new(idir) = sqrt(acell_new(idir))
512 [ + + ]: 26 : rprim(:,idir) = rprimd_new(:,idir)/acell_new(idir)
513 : : end do
514 : :
515 : : end if ! relaxstr == 1
516 : :
517 : : ! Write out the results
518 : :
519 [ + - ]: 2 : if (iwrite) then
520 : 2 : write(iout,*)
521 [ + - + + ]: 162 : write(iout,'(a,80a,a)') ch10,('=',ii=1,80),ch10
522 : 2 : write(iout,*)
523 : 2 : write(iout,*)'Relaxation of the geometry at fixed polarization:'
524 : 2 : write(iout,*)
525 : 2 : write(iout,'(a,3(2x,f16.9))')' Lambda = ',(lambda(idir),idir = 1, 3)
526 : 2 : write(iout,'(a,e16.9)')' Value of the energy functional E_1 = ',e1
527 : 2 : write(iout,*)
528 : 2 : write(iout,*)'Difference between actual value of the Polarization (C/m^2)'
529 : 2 : write(iout,*)'and the target value:'
530 : : end if
531 [ + + ]: 8 : diffpol(:) = (ptot_cart(:) - targetpol(:))*e_Cb/((Bohr_Ang*1.0d-10)**2)
532 [ + - ]: 2 : if (iwrite) write(iout,'(3(3x,f16.9))')(diffpol(idir),idir = 1, 3)
533 : :
534 [ + - ]: 2 : if (relaxat == 1) then
535 : : ! Compute the forces induced on the atoms by the electric field
536 : : ! The strength of the field is determined by lambda
537 [ + + + + ]: 30 : felfd(:,:) = zero
538 [ + + ]: 9 : do iatom = 1, natom
539 [ + + ]: 30 : do idir = 1, 3
540 [ + + ]: 91 : do jdir = 1, 3
541 : 63 : index = idir +3*((iatom - 1) + mpert*((jdir - 1) + 3*(natom + 1)))
542 : 84 : felfd(idir,iatom) = felfd(idir,iatom) - lambda(jdir)*blkval(1,index)/ucvol
543 : : end do
544 : : end do
545 : : end do
546 : :
547 : : ! Compute remaining forces and write them out
548 : :
549 [ + + + + ]: 30 : fdiff(:,:) = fcart(:,:) - felfd(:,:)
550 [ + - ]: 2 : if (iwrite) then
551 : 2 : write(iout,*)
552 : 2 : write(iout,*)'Difference between the Hellmann-Feynman forces'
553 : 2 : write(iout,*)'and the forces induced by the electric field'
554 : 2 : write(iout,*)'(cartesian coordinates, hartree/bohr)'
555 : : end if
556 : 2 : fmax = zero
557 [ + + ]: 9 : do iatom = 1, natom
558 [ + - ]: 7 : if (iwrite) write(iout,'(3(3x,es16.9))')(fdiff(idir,iatom),idir = 1, 3)
559 [ + + ]: 30 : do idir = 1, 3
560 [ + + ]: 28 : if (abs(fdiff(idir,iatom)) > fmax) fmax = abs(fdiff(idir,iatom))
561 : : end do
562 : : end do
563 : :
564 [ + - ]: 2 : if (iwrite) then
565 : 2 : write(iout,'(a,3x,es16.9)')' fmax = ',fmax
566 : 2 : write(iout,*)
567 : 2 : write(iout,*)'Change of cartesian coordinates (delta_xcart):'
568 [ + + ]: 9 : do iatom = 1, natom
569 : 9 : write(iout,'(5x,i3,3(2x,f16.9))')iatom,(delta_xcart(idir,iatom),idir = 1, 3)
570 : : end do
571 : 2 : write(iout,*)
572 : 2 : write(iout,*)'New cartesian coordinates (xcart_new):'
573 : 2 : write(iout,*)' xcart'
574 [ + + ]: 9 : do iatom = 1, natom
575 : 9 : write(iout,'(3(3x,d22.14))')(xcart_new(idir,iatom),idir = 1, 3)
576 : : end do
577 : 2 : write(iout,*)
578 : 2 : write(iout,*)'New reduced coordinates (xred_new):'
579 : 2 : write(iout,*)' xred'
580 [ + + ]: 9 : do iatom = 1, natom
581 : 9 : write(iout,'(3(3x,d22.14))')(xred_new(idir,iatom),idir = 1, 3)
582 : : end do
583 : : end if
584 : :
585 : : end if ! relaxat == 1
586 : :
587 [ + - ]: 2 : if (relaxstr == 1) then
588 : :
589 : : ! Compute the stresses induced by the electric field
590 : 2 : sigelfd(:) = zero
591 : 2 : istrain = 0
592 [ + + ]: 6 : do ipert = 1, 2
593 : 4 : jpert = natom + 2 + ipert
594 [ + + ]: 18 : do idir = 1, 3
595 : 12 : istrain = istrain + 1
596 [ + + ]: 48 : do jdir = 1, 3
597 : 36 : index = idir +3*((jpert - 1) + mpert*((jdir - 1) + 3*(natom + 1)))
598 : 48 : sigelfd(istrain) = sigelfd(istrain) + lambda(jdir)*blkval(1,index)
599 : : end do
600 : 16 : sigelfd(istrain) = sigelfd(istrain)/ucvol
601 : : end do
602 : : end do
603 : :
604 : : ! Compute the remaining stresses and write them out
605 [ + + ]: 14 : diffsig(:) = strten(:) - sigelfd(:)
606 : 2 : sigmax = zero
607 [ + + ]: 14 : do istrain = 1, 6
608 [ + + ]: 14 : if (abs(diffsig(istrain)) > sigmax) sigmax = abs(diffsig(istrain))
609 : : end do
610 [ + - ]: 2 : if (iwrite) then
611 : 2 : write(iout,*)
612 : 2 : write(iout,*)'Difference between the Hellmann-Feynman stresses'
613 : 2 : write(iout,*)'and the stresses induced by the electric field'
614 : 2 : write(iout,*)'(cartesian coordinates, hartree/bohr^3)'
615 : 2 : write(iout,'(2x,a,f16.9,5x,a,f16.9)')'diffsig(1) = ',diffsig(1),'diffsig(4) = ',diffsig(4)
616 : 2 : write(iout,'(2x,a,f16.9,5x,a,f16.9)')'diffsig(2) = ',diffsig(2),'diffsig(5) = ',diffsig(5)
617 : 2 : write(iout,'(2x,a,f16.9,5x,a,f16.9)')'diffsig(3) = ',diffsig(3),'diffsig(6) = ',diffsig(6)
618 : 2 : write(iout,'(a,3x,es16.9)')' sigmax = ',sigmax
619 : 2 : write(iout,*)
620 : 2 : write(iout,*)'Induced strain (delta_eta):'
621 : 2 : write(iout,'(2x,a,f16.9,5x,a,f16.9)')'delta_eta(1) = ',delta_eta(1),'delta_eta(4) = ',delta_eta(4)
622 : 2 : write(iout,'(2x,a,f16.9,5x,a,f16.9)')'delta_eta(2) = ',delta_eta(2),'delta_eta(5) = ',delta_eta(5)
623 : 2 : write(iout,'(2x,a,f16.9,5x,a,f16.9)')'delta_eta(3) = ',delta_eta(3),'delta_eta(6) = ',delta_eta(6)
624 : 2 : write(iout,*)
625 : 2 : write(iout,*)'New lattice constants (acell_new):'
626 : 2 : write(iout,*)' acell'
627 : 2 : write(iout,'(3(2x,d22.14))')(acell_new(idir),idir = 1, 3)
628 : 2 : write(iout,*)
629 : 2 : write(iout,*)'New primitive vectors (rprim_new):'
630 : 2 : write(iout,*)' rprim'
631 : 2 : write(iout,'(3(2x,d22.14))')(rprim(idir,1),idir = 1, 3)
632 : 2 : write(iout,'(3(2x,d22.14))')(rprim(idir,2),idir = 1, 3)
633 : 2 : write(iout,'(3(2x,d22.14))')(rprim(idir,3),idir = 1, 3)
634 : : end if
635 : : end if ! relaxstr /= 0
636 : :
637 : : end if ! (relaxat /= 0).or.(relaxstr /= 0)
638 : :
639 [ - + ]: 3 : ABI_FREE(cfac)
640 [ - + ]: 3 : ABI_FREE(fdiff)
641 [ - + ]: 3 : ABI_FREE(felfd)
642 [ - + ]: 3 : ABI_FREE(delta)
643 [ - + ]: 3 : ABI_FREE(fcart)
644 [ - + ]: 3 : ABI_FREE(fcmat)
645 [ - + ]: 3 : ABI_FREE(ifcmat)
646 [ - + ]: 3 : ABI_FREE(ipvt)
647 [ - + ]: 3 : ABI_FREE(rfpert)
648 [ - + ]: 3 : ABI_FREE(vec)
649 [ - + ]: 3 : ABI_FREE(zgwork)
650 [ - + ]: 3 : ABI_FREE(irelaxat)
651 : :
652 : 3 : end subroutine relaxpol
653 : : !!***
654 : :
655 : : end module m_relaxpol
656 : : !!***
|
