Line data Source code
1 : !!****m* ABINIT/m_ddb_flexo
2 : !! NAME
3 : !! m_ddb_flexo
4 : !!
5 : !! FUNCTION
6 : !! FIXME: add description.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2019-2026 ABINIT group (MR,MS)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! NOTES
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_ddb_flexo
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_profiling_abi
29 : use m_errors
30 :
31 : use m_fstrings, only : itoa,sjoin
32 : use m_ddb_hdr
33 : use m_ddb
34 : use m_crystal, only : crystal_t
35 : use m_dynmat, only : asria_corr,cart39
36 :
37 : implicit none
38 :
39 : private
40 :
41 : public :: ddb_flexo
42 : public :: ddb_phi1
43 : ! *************************************************************************
44 :
45 : contains
46 : !!***
47 :
48 : !!****f* ABINIT/m_ddb_flexo/ddb_flexo
49 : !! NAME
50 : !! ddb_flexo
51 : !!
52 : !! FUNCTION
53 : !! Get all the contributions to the flexoelectric tensor
54 : !!
55 : !! INPUTS
56 : !! asr= if /=0 acustic sume rule is imposed on the dynamical matrix
57 : !! d2asr(2,3,natom,3,natom)=ASR-correction
58 : !! ddb<type(ddb_type)>=2nd order derivative database.
59 : !! ddb_lw<type(ddb_type)>=Long wave 3rd order derivative database.
60 : !! ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files.
61 : !! Crystal<type(crystal_t)>=Crystal structure parameters
62 : !! filnamddb = name of the ddb file
63 : !! flexoflg= 1 -> Computes all contributions to FxE
64 : !! 2 -> Computes electronic (clamped ion) contribution to FxE
65 : !! 3 -> Computes mixed contribution to FxE
66 : !! 4 -> Computes lattice contribution to FxE
67 : !! prtvol= if > 1 print all individual quantities of the lattice contribution FxE
68 : !! zeff(3,3,natom)= Born Effective charges
69 : !!
70 : !! OUTPUT
71 : !!
72 : !! SIDE EFFECTS
73 : !!
74 : !! NOTES
75 : !!
76 : !! SOURCE
77 :
78 2 : subroutine ddb_flexo(asr,d2asr,ddb,ddb_lw,ddb_version,crystal,filnamddb,flexoflg,prtvol,zeff)
79 :
80 : !Arguments ------------------------------------
81 : !scalars
82 : integer,intent(in) :: ddb_version
83 : integer , intent(in) :: asr,flexoflg,prtvol
84 : class(ddb_type),intent(in) :: ddb,ddb_lw
85 : type(crystal_t),intent(in) :: crystal
86 : character(len=fnlen) :: filnamddb
87 : !arrays
88 : real(dp),intent(in) :: d2asr(2,3,ddb%natom,3,ddb%natom)
89 : real(dp),intent(in) :: zeff(3,3,ddb%natom)
90 :
91 : !Local variables-------------------------------
92 : integer :: elfd,iblok,ivar,jblok,kblok,lblok,lwsym,qvecd
93 : logical :: intstrn_only,iwrite
94 : character(len=500) :: msg
95 :
96 : !arrays
97 : integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/)
98 : integer :: rfelfd(4),rfphon(4),rfstrs(4)
99 : integer :: rfqvec(4)
100 : real(dp) :: qphnrm(3),qphon(3,3)
101 : real(dp) :: ciflexo(3,3,3,3)
102 4 : real(dp) :: intstrn(3,3,3,ddb%natom)
103 4 : real(dp) :: piezofr(3,ddb%natom,3,3)
104 : real(dp) :: lattflexo(3,3,3,3)
105 : real(dp) :: mixflexo(3,3,3,3)
106 4 : real(dp) :: pol1(3,3,3,ddb%natom)
107 4 : real(dp) :: psinvdm(3*ddb%natom,3*ddb%natom)
108 : real(dp) :: totflexo(3,3,3,3)
109 : character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/)
110 :
111 : ! *************************************************************************
112 :
113 : DBG_ENTER("COLL")
114 :
115 : ! First get the clamped-ion flexoelectric tensor
116 2 : ciflexo(:,:,:,:)=zero
117 2 : if (flexoflg==1.or.flexoflg==2) then
118 :
119 : rfphon(:)=0
120 2 : rfelfd(:)=0
121 2 : rfstrs(:)=0
122 2 : rfqvec(:)=0
123 :
124 : ! Look for the Gamma Block in the DDB
125 2 : qphon(:,:)=zero
126 8 : qphnrm(:)=one
127 2 : rfphon(:)=0
128 2 : rfelfd(1)=2
129 2 : rfstrs(2)=3
130 2 : rfqvec(3)=1
131 :
132 2 : write(msg, '(2a)' ) ch10," Extract the electronic flexoelectric coeficients from 3DTE"
133 2 : call wrtout(std_out,msg,'COLL')
134 2 : call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec)
135 :
136 2 : if (iblok == 0) then
137 0 : call wrtout(std_out, " ")
138 0 : call wrtout(std_out, "--- !WARNING")
139 0 : call wrtout(std_out, sjoin("- Cannot find clamped ion FxE tensor in DDB file:", filnamddb))
140 0 : call wrtout(std_out, " flexoflag=1 or 2 requires the DDB file to include the corresponding long wave 3rd derivatives")
141 : else
142 2 : call dtciflexo(ddb_lw%val(:,:,iblok),ddb_version,ddb%mpert,ddb%natom,ciflexo,crystal%ucvol)
143 : end if
144 :
145 : end if
146 :
147 : ! Then get the mixed contribution to the flexoelectric tensor
148 : !Activate the calculation of internal strain necessary for lattice mediated contribution.
149 2 : intstrn_only=.false.;if (flexoflg==4) intstrn_only=.true.
150 2 : mixflexo(:,:,:,:)=zero
151 2 : if (flexoflg==1.or.flexoflg==3.or.intstrn_only) then
152 :
153 : ! Extract the P^(1) tensor from the DDB
154 2 : if (.not.intstrn_only) then
155 2 : lwsym=0
156 2 : iblok = ddb_lw%get_quadrupoles(ddb_version,lwsym,BLKTYP_d3E_lw,pol1)
157 : end if
158 :
159 2 : rfphon(:)=0
160 2 : rfelfd(:)=0
161 2 : rfstrs(:)=0
162 2 : rfqvec(:)=0
163 :
164 : ! Look for the Gamma Block of the Phi^(1) tensor in the DDB
165 2 : qphon(:,:)=zero
166 8 : qphnrm(:)=one
167 2 : rfphon(1)=1
168 2 : rfphon(2)=1
169 2 : rfqvec(3)=1
170 :
171 2 : write(msg, '(2a)' ) ch10," Extract the Phi^(1) coeficients from 3DTE"
172 2 : call wrtout(std_out,msg,'COLL')
173 2 : call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec)
174 :
175 2 : if (iblok == 0) then
176 0 : call wrtout(std_out, " ")
177 0 : call wrtout(std_out, "--- !WARNING")
178 0 : call wrtout(std_out, sjoin("- Cannot find Phi^(1) tensor in DDB file:", filnamddb))
179 0 : call wrtout(std_out, " flexoflag=1 or 3 requires the DDB file to include the corresponding long wave 3rd derivatives")
180 : end if
181 :
182 : ! Look for th block that contains the forces
183 2 : qphon(:,:)=zero
184 8 : qphnrm(:)=one
185 2 : rfphon(:)=0
186 2 : rfphon(4)=1
187 :
188 2 : write(msg, '(2a)' ) ch10," Extract the forces from 1DTE"
189 2 : call wrtout(std_out,msg,'COLL')
190 2 : call ddb%get_block(kblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,4)
191 :
192 2 : if (kblok == 0) then
193 0 : call wrtout(std_out, " ")
194 0 : call wrtout(std_out, "--- !WARNING")
195 0 : call wrtout(std_out, sjoin("- Cannot find forces in DDB file:", filnamddb))
196 0 : call wrtout(std_out, " If there are nonzero residual atomic forces on the structure")
197 0 : call wrtout(std_out, " flexoflag=1 or 3 will produce an improper piezoelectric force response tensor")
198 0 : call wrtout(std_out, " and a wrong value for the mixed contribution to the flexoelectric tensor")
199 : end if
200 :
201 : ! Look for the Gamma Block of the dynamical matrix in the DDB
202 2 : qphon(:,:)=zero
203 8 : qphnrm(:)=one
204 2 : rfphon(:)=0
205 6 : rfphon(1:2)=1
206 :
207 2 : write(msg, '(2a)' ) ch10," Extract the Dynamical Matrix from 2DTE"
208 2 : call wrtout(std_out,msg,'COLL')
209 2 : call ddb%get_block(jblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,1)
210 :
211 2 : if (jblok == 0) then
212 0 : call wrtout(std_out, " ")
213 0 : call wrtout(std_out, "--- !WARNING")
214 0 : call wrtout(std_out, sjoin("- Cannot find Gamma point Dynamical Matrix in DDB file:", filnamddb))
215 0 : call wrtout(std_out, " flexoflag=1 or 3 requires the DDB file to include the corresponding 2nd derivatives")
216 : end if
217 :
218 2 : if (iblok/=0.and.jblok/=0) then
219 : call dtmixflexo(asr,d2asr,ddb%val(:,:,kblok),ddb%val(:,:,jblok),ddb_lw%val(:,:,iblok),ddb_version,crystal%gprimd,&
220 2 : & intstrn,intstrn_only,mixflexo,ddb%mpert,ddb%natom,piezofr,pol1,psinvdm,crystal%rprimd,crystal%ucvol)
221 : end if
222 :
223 : end if
224 :
225 : ! Finally get the lattice mediated contribution to the flexoelectric tensor
226 2 : lattflexo(:,:,:,:)=zero
227 2 : if (flexoflg==1.or.flexoflg==4) then
228 :
229 2 : rfphon(:)=0
230 2 : rfelfd(:)=0
231 2 : rfstrs(:)=0
232 2 : rfqvec(:)=0
233 :
234 : ! Look for the Gamma Block of the Phi^(1) tensor in the DDB
235 2 : qphon(:,:)=zero
236 8 : qphnrm(:)=one
237 2 : rfphon(1)=1
238 2 : rfphon(2)=1
239 2 : rfqvec(3)=1
240 :
241 2 : write(msg, '(2a)' ) ch10," Extract the Phi^(1) coeficients from 3DTE"
242 2 : call wrtout(std_out,msg,'COLL')
243 2 : call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec)
244 :
245 2 : if (iblok == 0) then
246 0 : call wrtout(std_out, " ")
247 0 : call wrtout(std_out, "--- !WARNING")
248 0 : call wrtout(std_out, sjoin("- Cannot find Phi^(1) tensor in DDB file:", filnamddb))
249 0 : call wrtout(std_out, " flexoflag=1 or 4 requires the DDB file to include the corresponding long wave 3rd derivatives")
250 : end if
251 :
252 : ! Look for the Gamma Block of the flexoelectric force response tensor in the DDB
253 2 : qphon(:,:)=zero
254 8 : qphnrm(:)=one
255 2 : rfphon(:)=0
256 2 : rfphon(1)=1
257 2 : rfstrs(2)=3
258 : rfqvec(3)=1
259 :
260 2 : write(msg, '(2a)' ) ch10," Extract the FxE force response coeficients from 3DTE"
261 2 : call wrtout(std_out,msg,'COLL')
262 2 : call ddb_lw%get_block(jblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec)
263 :
264 2 : if (jblok == 0) then
265 0 : call wrtout(std_out, " ")
266 0 : call wrtout(std_out, "--- !WARNING")
267 0 : call wrtout(std_out, sjoin("- Cannot find FxE force response tensor in DDB file:", filnamddb))
268 0 : call wrtout(std_out, " flexoflag=1 or 4 requires the DDB file to include the corresponding long wave 3rd derivatives")
269 : end if
270 :
271 : ! Look for the stress tensor in the DDB
272 2 : qphon(:,:)=zero
273 8 : qphnrm(:)=one
274 2 : rfphon(:)=0
275 2 : rfstrs(:)=0
276 2 : rfqvec(:)=0
277 2 : rfstrs(4)=3
278 :
279 2 : write(msg, '(2a)' ) ch10," Extract the stress tensor from 1DTE"
280 2 : call wrtout(std_out,msg,'COLL')
281 2 : call ddb%get_block(lblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,4)
282 :
283 2 : if (lblok == 0) then
284 0 : call wrtout(std_out, " ")
285 0 : call wrtout(std_out, "--- !WARNING")
286 0 : call wrtout(std_out, sjoin("- Cannot find stress tensor in DDB file:", filnamddb))
287 0 : call wrtout(std_out, " flexoflag=1 or 4 requires the DDB file to include the corresponding 2nd derivatives")
288 0 : call wrtout(std_out, " to compute the Lagrange Elastic Tensor ")
289 : end if
290 :
291 2 : if (iblok/=0.and.jblok/=0) then
292 : call dtlattflexo(ddb%amu,ddb%val(:,:,lblok),ddb_lw%val(:,:,jblok),ddb_lw%val(:,:,iblok),ddb_version,&
293 2 : & intstrn,lattflexo,ddb%mpert,ddb%natom,crystal%ntypat,piezofr,prtvol,psinvdm,crystal%typat,crystal%ucvol,zeff)
294 : end if
295 : end if
296 :
297 : !Merge the three contributions and print the total FxE tensor
298 242 : totflexo(:,:,:,:)=ciflexo(:,:,:,:)+mixflexo(:,:,:,:)+lattflexo(:,:,:,:)
299 :
300 2 : iwrite = ab_out > 0
301 2 : if (iwrite) then
302 2 : write(msg,'(3a)')ch10,' TOTAL flexoelectric tensor (units= nC/m) ',ch10
303 6 : call wrtout([ab_out,std_out],msg,'COLL')
304 2 : write(msg,*)' xx yy zz yz xz xy'
305 6 : call wrtout([ab_out,std_out],msg,'COLL')
306 14 : do ivar=1,6
307 12 : elfd=alpha(ivar)
308 12 : qvecd=beta(ivar)
309 12 : write(msg,'(3x,a2,6f12.6)') voigt(ivar),totflexo(elfd,qvecd,1,1),totflexo(elfd,qvecd,2,2),totflexo(elfd,qvecd,3,3),&
310 24 : totflexo(elfd,qvecd,2,3),totflexo(elfd,qvecd,1,3),totflexo(elfd,qvecd,1,2)
311 38 : call wrtout([ab_out,std_out],msg,'COLL')
312 : end do
313 8 : do ivar=4,6
314 6 : elfd=beta(ivar)
315 6 : qvecd=alpha(ivar)
316 6 : write(msg,'(3x,a2,6f12.6)') voigt(ivar+3),totflexo(elfd,qvecd,1,1),totflexo(elfd,qvecd,2,2),totflexo(elfd,qvecd,3,3),&
317 12 : totflexo(elfd,qvecd,2,3),totflexo(elfd,qvecd,1,3),totflexo(elfd,qvecd,1,2)
318 20 : call wrtout([ab_out,std_out],msg,'COLL')
319 : end do
320 : end if
321 :
322 :
323 : DBG_EXIT("COLL")
324 :
325 2 : end subroutine ddb_flexo
326 : !!***
327 :
328 : !!****f* m_ddb/dtciflexo
329 : !! NAME
330 : !! dtciflexo
331 : !!
332 : !! FUNCTION
333 : !! Reads the Clamped Ion Flexoelectric Tensor
334 : !! in the Gamma Block coming from the Derivative Data Base
335 : !! (long wave third-order derivatives).
336 : !!
337 : !! INPUTS
338 : !! blkval(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies
339 : !! mpert =maximum number of ipert
340 : !! natom= number of atoms in unit cell
341 : !! ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files.
342 : !! ucvol= unit cell volume
343 : !!
344 : !! OUTPUT
345 : !! ciflexo(3,3,3,3) = type-II Clamped Ion Flexoelectric Tensor
346 : !!
347 : !! SOURCE
348 :
349 2 : subroutine dtciflexo(blkval,ddb_version,mpert,natom,ciflexo,ucvol)
350 :
351 : !Arguments -------------------------------
352 : !scalars
353 : integer,intent(in) :: ddb_version,mpert,natom
354 : real(dp),intent(in) :: ucvol
355 : !arrays
356 : real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert)
357 : real(dp),intent(out) :: ciflexo(3,3,3,3)
358 :
359 : !Local variables -------------------------
360 : !scalars
361 : integer,parameter :: cvrsio8=20100401
362 : integer :: elfd,istrs,ivarA,ri,strsd,strsd1,strsd2,strst,qvecd
363 : logical :: iwrite
364 : real(dp) :: fac
365 : character(len=500) :: msg
366 : real(dp),parameter :: confac=e_Cb/Bohr_meter*1.d9
367 : !arrays
368 : integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/)
369 4 : real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert)
370 : character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/)
371 :
372 : ! *********************************************************************
373 :
374 : DBG_ENTER("COLL")
375 :
376 375606 : d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/))
377 375606 : d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/))
378 :
379 : !Define the factors to apply if DDB file has been created with the old version of
380 : !the longwave driver.
381 2 : if (ddb_version <= cvrsio8) then
382 1 : fac=-two/ucvol
383 1 : ri=2
384 : else
385 1 : fac=one/ucvol
386 1 : ri=1
387 : end if
388 :
389 : !Extraction of the clamped-ion flexoelectric coeficients
390 8 : do qvecd=1,3
391 44 : do istrs=1,6
392 36 : strsd1=alpha(istrs)
393 36 : strsd2=beta(istrs)
394 36 : strst=natom+3; if (istrs>3) strst=natom+4
395 18 : strsd=istrs; if (istrs>3) strsd=istrs-3
396 150 : do elfd=1,3
397 108 : ciflexo(elfd,qvecd,strsd1,strsd2)=fac*d3cart(ri,elfd,natom+2,strsd,strst,qvecd,natom+8)*confac
398 144 : if (istrs>3) ciflexo(elfd,qvecd,strsd2,strsd1)=ciflexo(elfd,qvecd,strsd1,strsd2)
399 : end do
400 : end do
401 : end do
402 :
403 : !Print results
404 2 : iwrite = ab_out > 0
405 2 : if (iwrite) then
406 2 : write(msg,'(3a)')ch10,' Type-II electronic (clamped ion) flexoelectric tensor (units= nC/m) ',ch10
407 6 : call wrtout([ab_out,std_out],msg,'COLL')
408 2 : write(msg,*)' xx yy zz yz xz xy'
409 6 : call wrtout([ab_out,std_out],msg,'COLL')
410 14 : do ivarA=1,6
411 12 : elfd=alpha(ivarA)
412 12 : qvecd=beta(ivarA)
413 12 : write(msg,'(3x,a2,6f12.6)') voigt(ivarA), ciflexo(elfd,qvecd,1,1),ciflexo(elfd,qvecd,2,2), &
414 12 : ciflexo(elfd,qvecd,3,3), ciflexo(elfd,qvecd,2,3), &
415 24 : ciflexo(elfd,qvecd,1,3),ciflexo(elfd,qvecd,1,2)
416 38 : call wrtout([ab_out,std_out],msg,'COLL')
417 : end do
418 8 : do ivarA=4,6
419 6 : elfd=beta(ivarA)
420 6 : qvecd=alpha(ivarA)
421 6 : write(msg,'(3x,a2,6f12.6)') voigt(ivarA+3), ciflexo(elfd,qvecd,1,1),ciflexo(elfd,qvecd,2,2), &
422 6 : ciflexo(elfd,qvecd,3,3), ciflexo(elfd,qvecd,2,3), &
423 12 : ciflexo(elfd,qvecd,1,3),ciflexo(elfd,qvecd,1,2)
424 20 : call wrtout([ab_out,std_out],msg,'COLL')
425 : end do
426 : end if
427 :
428 : DBG_EXIT("COLL")
429 :
430 2 : end subroutine dtciflexo
431 : !!***
432 :
433 : !!****f* m_ddb/dtmixflexo
434 : !! NAME
435 : !! dtmixflexo
436 : !!
437 : !! FUNCTION
438 : !! Reads the P^(1) and Phi^(1) tensors
439 : !! in the Gamma Block coming from the Derivative Data Base
440 : !! (long wave third-order derivatives). And computes the mixed
441 : !! contribution to the flexoelectric tensor.
442 : !!
443 : !! INPUTS
444 : !! asr= if /=0 acustic sume rule is imposed on the dynamical matrix
445 : !! d2asr(2,3,natom,3,natom)=ASR-correction
446 : !! blkval1d(2,3,mpert,3,mpert)= 1st derivative wrt atom displacements (at least)
447 : !! blkval2d(2,3,mpert,3,mpert)= 2nd derivatives wrt two atom displacements (at least)
448 : !! blkval(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies for Phi^(1) tensor
449 : !! ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files.
450 : !! gprimd(3,3)= basis vectors in the reciprocal space
451 : !! intstrn_only= activates only the calculation of the internal strain tensor
452 : !! mpert =maximum number of ipert
453 : !! natom= number of atoms in unit cell
454 : !! piezofr(3,natom,3,3)= piezoelectric force response tensor
455 : !! pol1(3,3,3,natom)= tensor with the polarization induced by an atomic displacement (P^(1))
456 : !! rprimd(3,3)= basis vectors in the real space
457 : !! ucvol= unit cell volume
458 : !!
459 : !! OUTPUT
460 : !! mixflexo(3,3,3,3) = type-II mixed contribution to the Flexoelectric Tensor
461 : !! intstrn(3,3,3,natom) = relaxed-ion internal strain tensor
462 : !! psinvdm(3*natom,3*natom) = pseudo inverse of dynamical matrix
463 : !!
464 : !! SOURCE
465 :
466 2 : subroutine dtmixflexo(asr,d2asr,blkval1d,blkval2d,blkval,ddb_version,gprimd,intstrn,intstrn_only, &
467 2 : & mixflexo,mpert,natom,piezofr,pol1,psinvdm,rprimd,ucvol)
468 :
469 : !Arguments -------------------------------
470 : !scalars
471 : integer,intent(in) :: asr,ddb_version,mpert,natom
472 : real(dp),intent(in) :: ucvol
473 : logical,intent(in) :: intstrn_only
474 : !arrays
475 : real(dp),intent(in) :: d2asr(2,3,mpert,3,mpert)
476 : real(dp),intent(in) :: blkval1d(2,3,mpert,3,mpert)
477 : real(dp),intent(in) :: blkval2d(2,3,mpert,3,mpert)
478 : real(dp),intent(in) :: blkval(2,3*mpert*3*mpert*3*mpert)
479 : real(dp),intent(in) :: gprimd(3,3)
480 : real(dp),intent(out) :: intstrn(3,3,3,natom)
481 : real(dp),intent(out) :: piezofr(3,natom,3,3)
482 : real(dp),intent(inout) :: pol1(3,3,3,natom)
483 : real(dp),intent(out) :: psinvdm(3*natom,3*natom)
484 : real(dp),intent(in) :: rprimd(3,3)
485 : real(dp),intent(out) :: mixflexo(3,3,3,3)
486 :
487 : !Local variables -------------------------
488 : !scalars
489 : integer,parameter :: cvrsio8=20100401
490 : integer :: elfd,iat,iatd,ivar,jat,jatd,jvar,katd,qvecd,qvecd2
491 : logical :: iwrite
492 : real(dp),parameter :: confac=e_Cb/Bohr_meter*1.d9
493 : real(dp) :: fac
494 : character(len=500) :: msg
495 : !arrays
496 : integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/)
497 4 : real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert)
498 4 : real(dp) :: redforces(3,natom),forces(3,natom)
499 4 : real(dp) :: phi1(3,natom,3,natom,3)
500 : integer :: flg1(3),flg2(3)
501 : real(dp) :: vec1(3),vec2(3)
502 : character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/)
503 :
504 : ! *********************************************************************
505 :
506 : DBG_ENTER("COLL")
507 :
508 375606 : d3cart(1,:,:,:,:,:,:) = reshape(blkval(1,:),shape = (/3,mpert,3,mpert,3,mpert/))
509 375606 : d3cart(2,:,:,:,:,:,:) = reshape(blkval(2,:),shape = (/3,mpert,3,mpert,3,mpert/))
510 :
511 : !P^(1) lacks the 1/ucvol factor
512 242 : pol1=pol1/ucvol
513 :
514 : !Define the factors to apply if DDB file has been created with the old version of
515 : !the longwave driver.
516 2 : if (ddb_version <= cvrsio8) then
517 : fac=-two
518 : else
519 1 : fac=-one
520 : end if
521 :
522 : !Extraction of Phi^(1) tensor
523 : phi1(:,:,:,:,:) = fac*RESHAPE(d3cart(2,1:3,1:natom,1:3,1:natom,1:3,natom+8), &
524 812 : & SHAPE=[3,natom,3,natom,3])
525 :
526 : !Extraction of the forces and conversion to cartesian coordinates
527 : !to acount for the improper contribution
528 26 : redforces(:,:)=-blkval1d(1,1:3,1:natom,1,1)
529 :
530 26 : forces(:,:)=redforces(:,:)
531 8 : flg1(:)=1
532 8 : do iat=1,natom
533 24 : vec1(:)=forces(:,iat)
534 6 : call cart39(flg1,flg2,gprimd,iat,natom,rprimd,vec1,vec2)
535 26 : forces(:,iat)=vec2(:)
536 : end do
537 :
538 : !Calculate the piezoelectric force-response tensor including the improper contribution
539 782 : piezofr(:,:,:,:) = SUM(phi1(1:3,1:natom,1:3,1:natom,1:3), DIM=4)
540 8 : do iatd=1,3
541 80 : piezofr(iatd,:,:,iatd) = piezofr(iatd,:,:,iatd) + TRANSPOSE(forces(:,:))
542 : end do
543 :
544 : !Calculate the ion-relaxed internal strain tensor
545 : !First we need to obtain the pseudo-inverse of the dynamical matrix
546 2 : call dm_psinv(asr,blkval2d,d2asr,ab_out,psinvdm,mpert,natom)
547 :
548 : !Perfom the product with the piezo force-response
549 242 : intstrn(:,:,:,:)=zero
550 8 : do qvecd=1,3
551 26 : do katd=1,3
552 78 : do iatd=1,3
553 234 : do iat=1,natom
554 162 : ivar=(iat-1)*3+iatd
555 702 : do jatd=1,3
556 2268 : do jat=1,natom
557 1620 : jvar=(jat-1)*3+jatd
558 :
559 : intstrn(qvecd,katd,iatd,iat)= intstrn(qvecd,katd,iatd,iat) + &
560 2106 : psinvdm(ivar,jvar)*piezofr(jatd,jat,katd,qvecd)
561 :
562 : end do
563 : end do
564 : end do
565 : end do
566 : end do
567 : end do
568 :
569 2 : if (.not.intstrn_only) then
570 : !Finally calculate the mixed contribution to the FxE tensor
571 2 : mixflexo(:,:,:,:)=zero
572 8 : do elfd=1,3
573 26 : do qvecd=1,3
574 78 : do katd=1,3
575 234 : do qvecd2=1,3
576 702 : do iatd=1,3
577 2106 : do iat=1,natom
578 :
579 : mixflexo(elfd,qvecd,katd,qvecd2)=mixflexo(elfd,qvecd,katd,qvecd2) - &
580 1944 : pol1(elfd,qvecd,iatd,iat)*intstrn(qvecd2,katd,iatd,iat)*confac
581 :
582 : end do
583 : end do
584 : end do
585 : end do
586 : end do
587 : end do
588 : end if
589 :
590 : !Print results
591 2 : iwrite = ab_out > 0
592 2 : if (iwrite) then
593 2 : write(msg,'(3a)')ch10,' Force-response internal strain tensor from long-wave magnitudes (units: Hartree/Bohr)',ch10
594 6 : call wrtout([ab_out,std_out],msg,'COLL')
595 2 : write(msg,*)' atom dir xx yy zz yz xz xy'
596 6 : call wrtout([ab_out,std_out],msg,'COLL')
597 8 : do iat=1,natom
598 6 : write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'x', piezofr(1,iat,1,1),piezofr(1,iat,2,2),piezofr(1,iat,3,3),&
599 12 : piezofr(1,iat,2,3),piezofr(1,iat,1,3),piezofr(1,iat,1,2)
600 18 : call wrtout([ab_out,std_out],msg,'COLL')
601 6 : write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'y', piezofr(2,iat,1,1),piezofr(2,iat,2,2),piezofr(2,iat,3,3),&
602 12 : piezofr(2,iat,2,3),piezofr(2,iat,1,3),piezofr(2,iat,1,2)
603 18 : call wrtout([ab_out,std_out],msg,'COLL')
604 6 : write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'z', piezofr(3,iat,1,1),piezofr(3,iat,2,2),piezofr(3,iat,3,3),&
605 12 : piezofr(3,iat,2,3),piezofr(3,iat,1,3),piezofr(3,iat,1,2)
606 20 : call wrtout([ab_out,std_out],msg,'COLL')
607 : end do
608 :
609 2 : write(msg,'(3a)')ch10,' Displacement-response internal strain tensor from long-wave magnitudes (units: Bohr)',ch10
610 6 : call wrtout([ab_out,std_out],msg,'COLL')
611 2 : write(msg,*)' atom dir xx yy zz yz xz xy'
612 6 : call wrtout([ab_out,std_out],msg,'COLL')
613 8 : do iat=1,natom
614 6 : write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'x', intstrn(1,1,1,iat),intstrn(2,2,1,iat),intstrn(3,3,1,iat),&
615 12 : intstrn(2,3,1,iat),intstrn(1,3,1,iat),intstrn(1,2,1,iat)
616 18 : call wrtout([ab_out,std_out],msg,'COLL')
617 6 : write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'y', intstrn(1,1,2,iat),intstrn(2,2,2,iat),intstrn(3,3,2,iat),&
618 12 : intstrn(2,3,2,iat),intstrn(1,3,2,iat),intstrn(1,2,2,iat)
619 18 : call wrtout([ab_out,std_out],msg,'COLL')
620 6 : write(msg,'(2x,i3,3x,a3,2x,6f12.6)') iat, 'z', intstrn(1,1,3,iat),intstrn(2,2,3,iat),intstrn(3,3,3,iat),&
621 12 : intstrn(2,3,3,iat),intstrn(1,3,3,iat),intstrn(1,2,3,iat)
622 20 : call wrtout([ab_out,std_out],msg,'COLL')
623 : end do
624 :
625 2 : if (.not.intstrn_only) then
626 2 : write(msg,'(3a)')ch10,' Type-II mixed contribution to flexoelectric tensor (units: nC/m)',ch10
627 6 : call wrtout([ab_out,std_out],msg,'COLL')
628 2 : write(msg,*)' xx yy zz yz xz xy'
629 6 : call wrtout([ab_out,std_out],msg,'COLL')
630 14 : do ivar=1,6
631 12 : elfd=alpha(ivar)
632 12 : qvecd=beta(ivar)
633 12 : write(msg,'(3x,a2,6f12.6)') voigt(ivar),mixflexo(elfd,qvecd,1,1),mixflexo(elfd,qvecd,2,2),mixflexo(elfd,qvecd,3,3),&
634 24 : mixflexo(elfd,qvecd,2,3),mixflexo(elfd,qvecd,1,3),mixflexo(elfd,qvecd,1,2)
635 38 : call wrtout([ab_out,std_out],msg,'COLL')
636 : end do
637 8 : do ivar=4,6
638 6 : elfd=beta(ivar)
639 6 : qvecd=alpha(ivar)
640 6 : write(msg,'(3x,a2,6f12.6)') voigt(ivar+3),mixflexo(elfd,qvecd,1,1),mixflexo(elfd,qvecd,2,2),mixflexo(elfd,qvecd,3,3),&
641 12 : mixflexo(elfd,qvecd,2,3),mixflexo(elfd,qvecd,1,3),mixflexo(elfd,qvecd,1,2)
642 20 : call wrtout([ab_out,std_out],msg,'COLL')
643 : end do
644 : end if
645 : end if
646 :
647 : DBG_EXIT("COLL")
648 :
649 2 : end subroutine dtmixflexo
650 : !!***
651 :
652 : !!****f* m_ddb/dtlattflexo
653 : !! NAME
654 : !! dtlattflexo
655 : !!
656 : !! FUNCTION
657 : !! Reads Phi^(1), flexoelectric force response and internal strain tensors
658 : !! in the Gamma Block coming from the Derivative Data Base
659 : !! (long wave third-order derivatives). And computes the lattice mediated
660 : !! contribution to the flexoelectric tensor.
661 : !! It also computes and writes the Lagrangian Elastic tensor.
662 : !!
663 : !! INPUTS
664 : !! amu(ntypat)=mass each atom type in the unit cell
665 : !! blkval1d(2,3,mpert,3,mpert)= 1st derivative wrt stress (at least)
666 : !! blkval2d(2,3,mpert,3,mpert)= 2nd derivatives wrt atom displacements and electric field (at least)
667 : !! blkvalA(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies for FxE force response tensor
668 : !! blkvalB(2,3*mpert*3*mpert*3*mpert)= matrix of third-order energies for Phi^(1) tensor
669 : !! ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files.
670 : !! intstrn(3,3,3,natom)= relaxed-ion internal strain tensor
671 : !! mpert= maximum number of ipert
672 : !! natom= number of atoms in unit cell
673 : !! piezofr(3,natom,3,3)= piezoelectric force response tensor (required to compute the Lagrange elastic tensor)
674 : !! prtvol= if >1 print all tensors entering the structure of lattflexo
675 : !! psinvdm(3*natom,3*natom) = pseudo inverse of dynamical matrix
676 : !! typat(natom)= Type of each atom in the unit cell
677 : !! ucvol= unit cell volume
678 : !!
679 : !! OUTPUT
680 : !! lattflexo(3,3,3,3) = type-II lattice contribution to the Flexoelectric Tensor
681 : !!
682 : !! SOURCE
683 :
684 2 : subroutine dtlattflexo(amu,blkval1d,blkvalA,blkvalB,ddb_version,intstrn,lattflexo,mpert,natom,&
685 2 : & ntypat,piezofr,prtvol,psinvdm,typat,ucvol,zeff)
686 :
687 : !Arguments -------------------------------
688 : !scalars
689 : integer,intent(in) :: ddb_version,mpert,natom,ntypat,prtvol
690 : real(dp),intent(in) :: ucvol
691 :
692 : !arrays
693 : integer,intent(in) :: typat(natom)
694 : real(dp),intent(in) :: amu(ntypat)
695 : real(dp),intent(in) :: blkval1d(2,3,mpert,3,mpert)
696 : real(dp),intent(in) :: blkvalA(2,3*mpert*3*mpert*3*mpert)
697 : real(dp),intent(in) :: blkvalB(2,3*mpert*3*mpert*3*mpert)
698 : real(dp),intent(in) :: intstrn(3,3,3,natom)
699 : real(dp),intent(in) :: piezofr(3,natom,3,3)
700 : real(dp),intent(in) :: psinvdm(3*natom,3*natom)
701 : real(dp),intent(in) :: zeff(3,3,natom)
702 : real(dp),intent(out) :: lattflexo(3,3,3,3)
703 :
704 : !Local variables -------------------------
705 : !scalars
706 : integer,parameter :: cvrsio8=20100401
707 : integer :: elfd,iat,iatd,istrs,ivar,jat,jatd,jvar,kat,katd,strsd
708 : integer :: strsd1,strsd2,strst,qvecd,qvecd2
709 : logical :: iwrite
710 : real(dp) :: fac,mtot
711 : real(dp),parameter :: confac=e_Cb/Bohr_meter*1.d9
712 : character(len=500) :: msg
713 : !arrays
714 : integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/)
715 : real(dp) :: frcelast_t2(3,3,3,3)
716 4 : real(dp) :: Csupkap(3,natom,3,3,3),Csupkapsum(3,3,3,3)
717 4 : real(dp) :: d3cart(2,3,mpert,3,mpert,3,mpert)
718 4 : real(dp) :: flexois(3,natom,3,3,3)
719 4 : real(dp) :: flexofr(3,natom,3,3,3)
720 4 : real(dp) :: hatCsupkap(3,natom,3,3,3)
721 : real(dp) :: lmcelast(3,3,3,3)
722 4 : real(dp) :: phi1(3,natom,3,natom,3)
723 : real(dp) :: ricelast_t2(3,3,3,3)
724 4 : real(dp) :: roundbkt_k(3,3,3,3,natom)
725 : real(dp) :: sqrbkt_t1(3,3,3,3)
726 : real(dp) :: stress(3,3)
727 : character(len=2) :: voigt(9)=(/'xx','yy','zz','yz','xz','xy','zy','zx','yx'/)
728 :
729 : ! MR: Kept for testing
730 : ! integer :: i,j,k,l
731 : ! real(dp) :: delik,deljk,delil,deljl
732 :
733 : ! *********************************************************************
734 :
735 : DBG_ENTER("COLL")
736 :
737 : !Extraction of the stress tensor
738 2 : stress(1,1)=blkval1d(1,1,natom+3,1,1)
739 2 : stress(2,2)=blkval1d(1,2,natom+3,1,1)
740 2 : stress(3,3)=blkval1d(1,3,natom+3,1,1)
741 2 : stress(2,3)=blkval1d(1,1,natom+4,1,1);stress(3,2)=stress(2,3)
742 2 : stress(1,3)=blkval1d(1,2,natom+4,1,1);stress(3,1)=stress(1,3)
743 2 : stress(1,2)=blkval1d(1,3,natom+4,1,1);stress(2,1)=stress(1,2)
744 :
745 : !Calculate the sublattice-dependent round bracket tensor of PRB 88,174106 (2013)
746 : !First we need to extract the Phi^(1) tensor
747 375606 : d3cart(1,:,:,:,:,:,:) = reshape(blkvalB(1,:),shape = (/3,mpert,3,mpert,3,mpert/))
748 375606 : d3cart(2,:,:,:,:,:,:) = reshape(blkvalB(2,:),shape = (/3,mpert,3,mpert,3,mpert/))
749 :
750 : !Define the factors to apply if DDB file has been created with the old version of
751 : !the longwave driver.
752 2 : if (ddb_version <= cvrsio8) then
753 : fac=-two
754 : else
755 1 : fac=-one
756 : end if
757 :
758 : phi1(:,:,:,:,:) = fac*RESHAPE(d3cart(2,1:3,1:natom,1:3,1:natom,1:3,natom+8), &
759 812 : & SHAPE=[3,natom,3,natom,3])
760 :
761 : !Now perform the multiplication with the internal strain
762 728 : roundbkt_k(:,:,:,:,:)=zero
763 8 : do iat=1,natom
764 26 : do iatd=1,3
765 78 : do qvecd2=1,3
766 234 : do jatd=1,3
767 702 : do qvecd=1,3
768 2106 : do katd=1,3
769 6804 : do kat=1,natom
770 : roundbkt_k(iatd,qvecd,jatd,qvecd2,iat)=roundbkt_k(iatd,qvecd,jatd,qvecd2,iat) &
771 6318 : + phi1(iatd,iat,katd,kat,qvecd)*intstrn(qvecd2,jatd,katd,kat)
772 : end do
773 : end do
774 : ! write(100,'(5i3,1x,f12.6)') iat, iatd,qvecd,jatd,qvecd2,roundbkt_k(iatd,qvecd,jatd,qvecd2,iat)
775 : end do
776 : end do
777 : end do
778 : end do
779 : end do
780 :
781 : !Calculate now the Lagrange elastic tensors
782 : !First we need to extract the clamped-ion flexoelectric force response tensor
783 375606 : d3cart(1,:,:,:,:,:,:) = reshape(blkvalA(1,:),shape = (/3,mpert,3,mpert,3,mpert/))
784 375606 : d3cart(2,:,:,:,:,:,:) = reshape(blkvalA(2,:),shape = (/3,mpert,3,mpert,3,mpert/))
785 :
786 : !Define the factors to apply if DDB file has been created with the old version of
787 : !the longwave driver.
788 2 : if (ddb_version <= cvrsio8) then
789 : fac=-two
790 : else
791 1 : fac=one
792 : end if
793 :
794 14 : do istrs=1,6
795 12 : strsd1=alpha(istrs)
796 12 : strsd2=beta(istrs)
797 12 : strst=natom+3; if (istrs>3) strst=natom+4
798 6 : strsd=istrs; if (istrs>3) strsd=istrs-3
799 50 : do qvecd=1,3
800 156 : do iat=1,natom
801 468 : do iatd=1,3
802 324 : flexofr(iatd,iat,qvecd,strsd1,strsd2)=fac*d3cart(1,iatd,iat,strsd,strst,qvecd,natom+8)
803 432 : if (istrs>3) flexofr(iatd,iat,qvecd,strsd2,strsd1)=flexofr(iatd,iat,qvecd,strsd1,strsd2)
804 : end do
805 : end do
806 : end do
807 : end do
808 :
809 : !Now compute the type-II frozen-ion elastic tensor (without stress corrected)
810 2 : fac=one/ucvol
811 728 : frcelast_t2(:,:,:,:) = fac*SUM(flexofr(1:3,1:natom,1:3,1:3,1:3), DIM=2)
812 :
813 : !Now convert to type-I to obtain the square bracketed tensor of Born and Huang
814 8 : do qvecd=1,3
815 26 : do strsd2=1,3
816 78 : do strsd1=1,3
817 : sqrbkt_t1(1:3,strsd1,strsd2,qvecd)=half*(frcelast_t2(1:3,qvecd,strsd1,strsd2) + &
818 234 : & frcelast_t2(1:3,strsd2,strsd1,qvecd))
819 : end do
820 : end do
821 : end do
822 :
823 : !Now correct the stress in the square bracketed tesnor tensor
824 8 : do qvecd=1,3
825 26 : do strsd2=1,3
826 78 : do strsd1=1,3
827 234 : do iatd=1,3
828 216 : if (iatd==strsd1) then
829 54 : sqrbkt_t1(iatd,strsd1,strsd2,qvecd)=sqrbkt_t1(iatd,strsd1,strsd2,qvecd) - stress(strsd2,qvecd)
830 : endif
831 : end do
832 : end do
833 : end do
834 : end do
835 :
836 : !Now convert back to type-II in order to obtain the frozen ion Lagrange elastic tensor.
837 2 : frcelast_t2(:,:,:,:)=zero
838 8 : do iatd=1,3
839 26 : do qvecd=1,3
840 78 : do strsd1=1,3
841 234 : do strsd2=1,3
842 : frcelast_t2(iatd,qvecd,strsd1,strsd2)=sqrbkt_t1(iatd,strsd1,qvecd,strsd2) + &
843 216 : & sqrbkt_t1(iatd,strsd2,strsd1,qvecd)-sqrbkt_t1(iatd,qvecd,strsd2,strsd1)
844 : end do
845 : end do
846 : end do
847 : end do
848 :
849 : !MR: kept for testing only. If uncommented the resulting elastic tensors must agree with HWRV's ones
850 : ! write(ab_out,*)
851 : ! write(ab_out,*) "Stress tensor"
852 : ! do i=1,3
853 : ! write(ab_out,*) stress(i,:)
854 : ! end do
855 : ! do i=1,3
856 : ! do j=1,3
857 : ! do k=1,3
858 : ! delik=0.d0; if(i==k) delik=1.d0
859 : ! deljk=0.d0; if(j==k) deljk=1.d0
860 : ! do l=1,3
861 : ! delil=0.d0; if(i==l) delil=1.d0
862 : ! deljl=0.d0; if(j==l) deljl=1.d0
863 :
864 : ! frcelast_t2(i,j,k,l)=frcelast_t2(i,j,k,l) + &
865 : ! & 0.5d0*( deljk*stress(i,l) + delik*stress(j,l) + delil*stress(j,k) &
866 : ! & + deljl*stress(i,k) )
867 :
868 : ! end do
869 : ! end do
870 : ! end do
871 : ! end do
872 :
873 : !Now compute the contribution to the elastic tensor due to ion relaxations
874 : !and sum with the clamped ion elastic tensor to obtain the relaxed ion one
875 2 : lmcelast(:,:,:,:)=zero
876 8 : do strsd2=1,3
877 26 : do strsd1=1,3
878 78 : do qvecd=1,3
879 234 : do jatd=1,3
880 702 : do iatd=1,3
881 2106 : do iat=1,natom
882 : lmcelast(jatd,qvecd,strsd1,strsd2)=lmcelast(jatd,qvecd,strsd1,strsd2) - &
883 1944 : intstrn(jatd,qvecd,iatd,iat)*piezofr(iatd,iat,strsd1,strsd2)*fac
884 : end do
885 : end do
886 : end do
887 : end do
888 : end do
889 : end do
890 242 : ricelast_t2(:,:,:,:) = frcelast_t2(:,:,:,:) + lmcelast(:,:,:,:)
891 :
892 : !In last place compute the lattice contribution to the FxE tensor
893 : !First obtain the C^{\kappa} tensor of Eq. 59 of PRB 88,174106 (2013)
894 : !and its sublattice summation
895 2 : Csupkapsum(:,:,:,:)=zero
896 8 : do strsd2=1,3
897 26 : do strsd1=1,3
898 78 : do qvecd=1,3
899 234 : do iatd=1,3
900 702 : do iat=1,natom
901 : Csupkap(iatd,iat,qvecd,strsd1,strsd2)=flexofr(iatd,iat,qvecd,strsd1,strsd2) + &
902 486 : & roundbkt_k(iatd,qvecd,strsd1,strsd2,iat)
903 : Csupkapsum(iatd,qvecd,strsd1,strsd2)=Csupkapsum(iatd,qvecd,strsd1,strsd2) + &
904 648 : & Csupkap(iatd,iat,qvecd,strsd1,strsd2)
905 : end do
906 : end do
907 : end do
908 : end do
909 : end do
910 :
911 : !Then separate the mass-dependent part
912 2 : mtot=zero
913 8 : do iat=1,natom
914 8 : mtot=mtot + amu(typat(iat))
915 : end do
916 :
917 8 : do strsd2=1,3
918 26 : do strsd1=1,3
919 78 : do qvecd=1,3
920 234 : do iatd=1,3
921 702 : do iat=1,natom
922 : hatCsupkap(iatd,iat,qvecd,strsd1,strsd2)=Csupkap(iatd,iat,qvecd,strsd1,strsd2) - &
923 648 : & amu(typat(iat))/mtot*Csupkapsum(iatd,qvecd,strsd1,strsd2)
924 : end do
925 : end do
926 : end do
927 : end do
928 : end do
929 :
930 : !Now compute the type-II flexoelectric internal strain tensor
931 728 : flexois(:,:,:,:,:)=zero
932 8 : do strsd2=1,3
933 26 : do strsd1=1,3
934 78 : do qvecd=1,3
935 234 : do iat=1,natom
936 702 : do iatd=1,3
937 486 : ivar=(iat-1)*3+iatd
938 2268 : do jat=1,natom
939 6966 : do jatd=1,3
940 4860 : jvar=(jat-1)*3+jatd
941 : flexois(iatd,iat,qvecd,strsd1,strsd2)=flexois(iatd,iat,qvecd,strsd1,strsd2) + &
942 6480 : & psinvdm(ivar,jvar)*hatCsupkap(jatd,jat,qvecd,strsd1,strsd2)
943 : end do
944 : end do
945 : end do
946 : end do
947 : end do
948 : end do
949 : end do
950 :
951 : !Finally multiply by the effective charges to obtain the FxE tensor
952 2 : lattflexo(:,:,:,:)=zero
953 8 : do strsd2=1,3
954 26 : do strsd1=1,3
955 78 : do qvecd=1,3
956 234 : do elfd=1,3
957 702 : do iat=1,natom
958 2106 : do iatd=1,3
959 : lattflexo(elfd,qvecd,strsd1,strsd2)=lattflexo(elfd,qvecd,strsd1,strsd2) + &
960 1944 : & zeff(elfd,iatd,iat)*flexois(iatd,iat,qvecd,strsd1,strsd2)/ucvol*confac
961 : end do
962 : end do
963 : end do
964 : end do
965 : end do
966 : end do
967 :
968 :
969 2 : iwrite = ab_out > 0
970 2 : if (iwrite) then
971 2 : write(msg,'(3a)')ch10,' Lagrange elastic tensor from long wave magnitudes (clamped ion) (units= 10^2 GPa) ',ch10
972 6 : call wrtout([ab_out,std_out],msg,'COLL')
973 2 : write(msg,*)' xx yy zz yz xz xy'
974 6 : call wrtout([ab_out,std_out],msg,'COLL')
975 242 : frcelast_t2(:,:,:,:)=frcelast_t2(:,:,:,:)*HaBohr3_GPa/100.00_dp
976 14 : do ivar=1,6
977 12 : iatd=alpha(ivar)
978 12 : qvecd=beta(ivar)
979 12 : write(msg,'(9f12.6)') frcelast_t2(iatd,qvecd,1,1), frcelast_t2(iatd,qvecd,2,2),frcelast_t2(iatd,qvecd,3,3), &
980 24 : frcelast_t2(iatd,qvecd,2,3), frcelast_t2(iatd,qvecd,1,3),frcelast_t2(iatd,qvecd,1,2)
981 :
982 38 : call wrtout([ab_out,std_out],msg,'COLL')
983 : end do
984 :
985 2 : write(msg,'(3a)')ch10,' Lagrange elastic tensor from long wave magnitudes (relaxed ion) (units= 10^2 GPa) ',ch10
986 6 : call wrtout([ab_out,std_out],msg,'COLL')
987 2 : write(msg,*)' xx yy zz yz xz xy'
988 6 : call wrtout([ab_out,std_out],msg,'COLL')
989 242 : ricelast_t2(:,:,:,:)=ricelast_t2(:,:,:,:)*HaBohr3_GPa/100.00_dp
990 14 : do ivar=1,6
991 12 : iatd=alpha(ivar)
992 12 : qvecd=beta(ivar)
993 12 : write(msg,'(9f12.6)') ricelast_t2(iatd,qvecd,1,1), ricelast_t2(iatd,qvecd,2,2),ricelast_t2(iatd,qvecd,3,3), &
994 24 : ricelast_t2(iatd,qvecd,2,3), ricelast_t2(iatd,qvecd,1,3),ricelast_t2(iatd,qvecd,1,2)
995 :
996 38 : call wrtout([ab_out,std_out],msg,'COLL')
997 : end do
998 2 : if (prtvol > 1) then
999 0 : write(msg,'(3a)')ch10,' (...)^kappa contribution to the flexoelectric force-response tensor (units: eV)',ch10
1000 0 : call wrtout([ab_out,std_out],msg,'COLL')
1001 0 : write(msg,*)' atom dir xx yy zz yz xz xy'
1002 0 : call wrtout([ab_out,std_out],msg,'COLL')
1003 0 : roundbkt_k(:,:,:,:,:)=roundbkt_k(:,:,:,:,:)*Ha_eV
1004 0 : do iat=1,natom
1005 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', roundbkt_k(1,1,1,1,iat), roundbkt_k(1,1,2,2,iat), &
1006 0 : & roundbkt_k(1,1,3,3,iat),roundbkt_k(1,1,2,3,iat),roundbkt_k(1,1,1,3,iat),roundbkt_k(1,1,1,2,iat)
1007 0 : call wrtout([ab_out,std_out],msg,'COLL')
1008 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', roundbkt_k(2,2,1,1,iat), roundbkt_k(2,2,2,2,iat), &
1009 0 : & roundbkt_k(2,2,3,3,iat),roundbkt_k(2,2,2,3,iat),roundbkt_k(2,2,1,3,iat),roundbkt_k(2,2,1,2,iat)
1010 0 : call wrtout([ab_out,std_out],msg,'COLL')
1011 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', roundbkt_k(3,3,1,1,iat), roundbkt_k(3,3,2,2,iat), &
1012 0 : & roundbkt_k(3,3,3,3,iat),roundbkt_k(3,3,2,3,iat),roundbkt_k(3,3,1,3,iat),roundbkt_k(3,3,1,2,iat)
1013 0 : call wrtout([ab_out,std_out],msg,'COLL')
1014 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', roundbkt_k(2,3,1,1,iat), roundbkt_k(2,3,2,2,iat), &
1015 0 : & roundbkt_k(2,3,3,3,iat),roundbkt_k(2,3,2,3,iat),roundbkt_k(2,3,1,3,iat),roundbkt_k(2,3,1,2,iat)
1016 0 : call wrtout([ab_out,std_out],msg,'COLL')
1017 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', roundbkt_k(1,3,1,1,iat), roundbkt_k(1,3,2,2,iat), &
1018 0 : & roundbkt_k(1,3,3,3,iat),roundbkt_k(1,3,2,3,iat),roundbkt_k(1,3,1,3,iat),roundbkt_k(1,3,1,2,iat)
1019 0 : call wrtout([ab_out,std_out],msg,'COLL')
1020 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', roundbkt_k(1,2,1,1,iat), roundbkt_k(1,2,2,2,iat), &
1021 0 : & roundbkt_k(1,2,3,3,iat),roundbkt_k(1,2,2,3,iat),roundbkt_k(1,2,1,3,iat),roundbkt_k(1,2,1,2,iat)
1022 0 : call wrtout([ab_out,std_out],msg,'COLL')
1023 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', roundbkt_k(3,2,1,1,iat), roundbkt_k(3,2,2,2,iat), &
1024 0 : & roundbkt_k(3,2,3,3,iat),roundbkt_k(3,2,2,3,iat),roundbkt_k(3,2,1,3,iat),roundbkt_k(3,2,1,2,iat)
1025 0 : call wrtout([ab_out,std_out],msg,'COLL')
1026 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', roundbkt_k(3,1,1,1,iat), roundbkt_k(3,1,2,2,iat), &
1027 0 : & roundbkt_k(3,1,3,3,iat),roundbkt_k(3,1,2,3,iat),roundbkt_k(3,1,1,3,iat),roundbkt_k(3,1,1,2,iat)
1028 0 : call wrtout([ab_out,std_out],msg,'COLL')
1029 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', roundbkt_k(2,1,1,1,iat), roundbkt_k(2,1,2,2,iat), &
1030 0 : & roundbkt_k(2,1,3,3,iat),roundbkt_k(2,1,2,3,iat),roundbkt_k(2,1,1,3,iat),roundbkt_k(2,1,1,2,iat)
1031 0 : call wrtout([ab_out,std_out],msg,'COLL')
1032 : end do
1033 :
1034 0 : write(msg,'(3a)')ch10,' [...]^kappa contribution to the flexoelectric force-response tensor (units: eV)',ch10
1035 0 : call wrtout([ab_out,std_out],msg,'COLL')
1036 0 : write(msg,*)' atom dir xx yy zz yz xz xy'
1037 0 : call wrtout([ab_out,std_out],msg,'COLL')
1038 0 : flexofr(:,:,:,:,:)=flexofr(:,:,:,:,:)*Ha_eV
1039 0 : do iat=1,natom
1040 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', flexofr(1,iat,1,1,1),flexofr(1,iat,1,2,2),flexofr(1,iat,1,3,3),&
1041 0 : & flexofr(1,iat,1,2,3),flexofr(1,iat,1,1,3),flexofr(1,iat,1,1,2)
1042 0 : call wrtout([ab_out,std_out],msg,'COLL')
1043 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', flexofr(2,iat,2,1,1),flexofr(2,iat,2,2,2),flexofr(2,iat,2,3,3),&
1044 0 : & flexofr(2,iat,2,2,3),flexofr(2,iat,2,1,3),flexofr(2,iat,2,1,2)
1045 0 : call wrtout([ab_out,std_out],msg,'COLL')
1046 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', flexofr(3,iat,3,1,1),flexofr(3,iat,3,2,2),flexofr(3,iat,3,3,3),&
1047 0 : & flexofr(3,iat,3,2,3),flexofr(3,iat,3,1,3),flexofr(3,iat,3,1,2)
1048 0 : call wrtout([ab_out,std_out],msg,'COLL')
1049 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', flexofr(2,iat,3,1,1),flexofr(2,iat,3,2,2),flexofr(2,iat,3,3,3),&
1050 0 : & flexofr(2,iat,3,2,3),flexofr(2,iat,3,1,3),flexofr(2,iat,3,1,2)
1051 0 : call wrtout([ab_out,std_out],msg,'COLL')
1052 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', flexofr(1,iat,3,1,1),flexofr(1,iat,3,2,2),flexofr(1,iat,3,3,3),&
1053 0 : & flexofr(1,iat,3,2,3),flexofr(1,iat,3,1,3),flexofr(1,iat,3,1,2)
1054 0 : call wrtout([ab_out,std_out],msg,'COLL')
1055 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', flexofr(1,iat,2,1,1),flexofr(1,iat,2,2,2),flexofr(1,iat,2,3,3),&
1056 0 : & flexofr(1,iat,2,2,3),flexofr(1,iat,2,1,3),flexofr(1,iat,2,1,2)
1057 0 : call wrtout([ab_out,std_out],msg,'COLL')
1058 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', flexofr(3,iat,2,1,1),flexofr(3,iat,2,2,2),flexofr(3,iat,2,3,3),&
1059 0 : & flexofr(3,iat,2,2,3),flexofr(3,iat,2,1,3),flexofr(3,iat,2,1,2)
1060 0 : call wrtout([ab_out,std_out],msg,'COLL')
1061 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', flexofr(3,iat,1,1,1),flexofr(3,iat,1,2,2),flexofr(3,iat,1,3,3),&
1062 0 : & flexofr(3,iat,1,2,3),flexofr(3,iat,1,1,3),flexofr(3,iat,1,1,2)
1063 0 : call wrtout([ab_out,std_out],msg,'COLL')
1064 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', flexofr(2,iat,1,1,1),flexofr(2,iat,1,2,2),flexofr(2,iat,1,3,3),&
1065 0 : & flexofr(2,iat,1,2,3),flexofr(2,iat,1,1,3),flexofr(2,iat,1,1,2)
1066 0 : call wrtout([ab_out,std_out],msg,'COLL')
1067 : end do
1068 : end if
1069 :
1070 2 : write(msg,'(3a)')ch10,' Flexoelectric force-response tensor (units: eV)',ch10
1071 6 : call wrtout([ab_out,std_out],msg,'COLL')
1072 2 : write(msg,*)' atom dir xx yy zz yz xz xy'
1073 6 : call wrtout([ab_out,std_out],msg,'COLL')
1074 728 : Csupkap(:,:,:,:,:)=Csupkap(:,:,:,:,:)*Ha_eV
1075 8 : do iat=1,natom
1076 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', Csupkap(1,iat,1,1,1),Csupkap(1,iat,1,2,2),Csupkap(1,iat,1,3,3),&
1077 12 : & Csupkap(1,iat,1,2,3),Csupkap(1,iat,1,1,3),Csupkap(1,iat,1,1,2)
1078 18 : call wrtout([ab_out,std_out],msg,'COLL')
1079 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', Csupkap(2,iat,2,1,1),Csupkap(2,iat,2,2,2),Csupkap(2,iat,2,3,3),&
1080 12 : & Csupkap(2,iat,2,2,3),Csupkap(2,iat,2,1,3),Csupkap(2,iat,2,1,2)
1081 18 : call wrtout([ab_out,std_out],msg,'COLL')
1082 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', Csupkap(3,iat,3,1,1),Csupkap(3,iat,3,2,2),Csupkap(3,iat,3,3,3),&
1083 12 : & Csupkap(3,iat,3,2,3),Csupkap(3,iat,3,1,3),Csupkap(3,iat,3,1,2)
1084 18 : call wrtout([ab_out,std_out],msg,'COLL')
1085 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', Csupkap(2,iat,3,1,1),Csupkap(2,iat,3,2,2),Csupkap(2,iat,3,3,3),&
1086 12 : & Csupkap(2,iat,3,2,3),Csupkap(2,iat,3,1,3),Csupkap(2,iat,3,1,2)
1087 18 : call wrtout([ab_out,std_out],msg,'COLL')
1088 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', Csupkap(1,iat,3,1,1),Csupkap(1,iat,3,2,2),Csupkap(1,iat,3,3,3),&
1089 12 : & Csupkap(1,iat,3,2,3),Csupkap(1,iat,3,1,3),Csupkap(1,iat,3,1,2)
1090 18 : call wrtout([ab_out,std_out],msg,'COLL')
1091 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', Csupkap(1,iat,2,1,1),Csupkap(1,iat,2,2,2),Csupkap(1,iat,2,3,3),&
1092 12 : & Csupkap(1,iat,2,2,3),Csupkap(1,iat,2,1,3),Csupkap(1,iat,2,1,2)
1093 18 : call wrtout([ab_out,std_out],msg,'COLL')
1094 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', Csupkap(3,iat,2,1,1),Csupkap(3,iat,2,2,2),Csupkap(3,iat,2,3,3),&
1095 12 : & Csupkap(3,iat,2,2,3),Csupkap(3,iat,2,1,3),Csupkap(3,iat,2,1,2)
1096 18 : call wrtout([ab_out,std_out],msg,'COLL')
1097 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', Csupkap(3,iat,1,1,1),Csupkap(3,iat,1,2,2),Csupkap(3,iat,1,3,3),&
1098 12 : & Csupkap(3,iat,1,2,3),Csupkap(3,iat,1,1,3),Csupkap(3,iat,1,1,2)
1099 18 : call wrtout([ab_out,std_out],msg,'COLL')
1100 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', Csupkap(2,iat,1,1,1),Csupkap(2,iat,1,2,2),Csupkap(2,iat,1,3,3),&
1101 12 : & Csupkap(2,iat,1,2,3),Csupkap(2,iat,1,1,3),Csupkap(2,iat,1,1,2)
1102 20 : call wrtout([ab_out,std_out],msg,'COLL')
1103 : end do
1104 2 : if (prtvol > 1) then
1105 0 : write(msg,'(3a)')ch10,' Flexoelectric force-response tensor minus mass dependent part (units: eV)',ch10
1106 0 : call wrtout([ab_out,std_out],msg,'COLL')
1107 0 : write(msg,*)' atom dir xx yy zz yz xz xy'
1108 0 : call wrtout([ab_out,std_out],msg,'COLL')
1109 0 : hatCsupkap(:,:,:,:,:)=hatCsupkap(:,:,:,:,:)*Ha_eV
1110 0 : do iat=1,natom
1111 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'xx',hatCsupkap(1,iat,1,1,1),hatCsupkap(1,iat,1,2,2),hatCsupkap(1,iat,1,3,3),&
1112 0 : & hatCsupkap(1,iat,1,2,3),hatCsupkap(1,iat,1,1,3),hatCsupkap(1,iat,1,1,2)
1113 0 : call wrtout([ab_out,std_out],msg,'COLL')
1114 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'yy',hatCsupkap(2,iat,2,1,1),hatCsupkap(2,iat,2,2,2),hatCsupkap(2,iat,2,3,3),&
1115 0 : & hatCsupkap(2,iat,2,2,3),hatCsupkap(2,iat,2,1,3),hatCsupkap(2,iat,2,1,2)
1116 0 : call wrtout([ab_out,std_out],msg,'COLL')
1117 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'zz',hatCsupkap(3,iat,3,1,1),hatCsupkap(3,iat,3,2,2),hatCsupkap(3,iat,3,3,3),&
1118 0 : & hatCsupkap(3,iat,3,2,3),hatCsupkap(3,iat,3,1,3),hatCsupkap(3,iat,3,1,2)
1119 0 : call wrtout([ab_out,std_out],msg,'COLL')
1120 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'yz',hatCsupkap(2,iat,3,1,1),hatCsupkap(2,iat,3,2,2),hatCsupkap(2,iat,3,3,3),&
1121 0 : & hatCsupkap(2,iat,3,2,3),hatCsupkap(2,iat,3,1,3),hatCsupkap(2,iat,3,1,2)
1122 0 : call wrtout([ab_out,std_out],msg,'COLL')
1123 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'xz',hatCsupkap(1,iat,3,1,1),hatCsupkap(1,iat,3,2,2),hatCsupkap(1,iat,3,3,3),&
1124 0 : & hatCsupkap(1,iat,3,2,3),hatCsupkap(1,iat,3,1,3),hatCsupkap(1,iat,3,1,2)
1125 0 : call wrtout([ab_out,std_out],msg,'COLL')
1126 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'xy',hatCsupkap(1,iat,2,1,1),hatCsupkap(1,iat,2,2,2),hatCsupkap(1,iat,2,3,3),&
1127 0 : & hatCsupkap(1,iat,2,2,3),hatCsupkap(1,iat,2,1,3),hatCsupkap(1,iat,2,1,2)
1128 0 : call wrtout([ab_out,std_out],msg,'COLL')
1129 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'zy',hatCsupkap(3,iat,2,1,1),hatCsupkap(3,iat,2,2,2),hatCsupkap(3,iat,2,3,3),&
1130 0 : & hatCsupkap(3,iat,2,2,3),hatCsupkap(3,iat,2,1,3),hatCsupkap(3,iat,2,1,2)
1131 0 : call wrtout([ab_out,std_out],msg,'COLL')
1132 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'zx',hatCsupkap(3,iat,1,1,1),hatCsupkap(3,iat,1,2,2),hatCsupkap(3,iat,1,3,3),&
1133 0 : & hatCsupkap(3,iat,1,2,3),hatCsupkap(3,iat,1,1,3),hatCsupkap(3,iat,1,1,2)
1134 0 : call wrtout([ab_out,std_out],msg,'COLL')
1135 0 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat,'yx',hatCsupkap(2,iat,1,1,1),hatCsupkap(2,iat,1,2,2),hatCsupkap(2,iat,1,3,3),&
1136 0 : & hatCsupkap(2,iat,1,2,3),hatCsupkap(2,iat,1,1,3),hatCsupkap(2,iat,1,1,2)
1137 0 : call wrtout([ab_out,std_out],msg,'COLL')
1138 : end do
1139 : end if
1140 :
1141 2 : write(msg,'(3a)')ch10,' Displacement-response flexoelectric internal strain tensor (units: Bohr^2)',ch10
1142 6 : call wrtout([ab_out,std_out],msg,'COLL')
1143 2 : write(msg,*)' atom dir xx yy zz yz xz xy'
1144 6 : call wrtout([ab_out,std_out],msg,'COLL')
1145 8 : do iat=1,natom
1146 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xx', flexois(1,iat,1,1,1),flexois(1,iat,1,2,2),flexois(1,iat,1,3,3),&
1147 12 : & flexois(1,iat,1,2,3),flexois(1,iat,1,1,3),flexois(1,iat,1,1,2)
1148 18 : call wrtout([ab_out,std_out],msg,'COLL')
1149 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yy', flexois(2,iat,2,1,1),flexois(2,iat,2,2,2),flexois(2,iat,2,3,3),&
1150 12 : & flexois(2,iat,2,2,3),flexois(2,iat,2,1,3),flexois(2,iat,2,1,2)
1151 18 : call wrtout([ab_out,std_out],msg,'COLL')
1152 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zz', flexois(3,iat,3,1,1),flexois(3,iat,3,2,2),flexois(3,iat,3,3,3),&
1153 12 : & flexois(3,iat,3,2,3),flexois(3,iat,3,1,3),flexois(3,iat,3,1,2)
1154 18 : call wrtout([ab_out,std_out],msg,'COLL')
1155 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yz', flexois(2,iat,3,1,1),flexois(2,iat,3,2,2),flexois(2,iat,3,3,3),&
1156 12 : & flexois(2,iat,3,2,3),flexois(2,iat,3,1,3),flexois(2,iat,3,1,2)
1157 18 : call wrtout([ab_out,std_out],msg,'COLL')
1158 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xz', flexois(1,iat,3,1,1),flexois(1,iat,3,2,2),flexois(1,iat,3,3,3),&
1159 12 : & flexois(1,iat,3,2,3),flexois(1,iat,3,1,3),flexois(1,iat,3,1,2)
1160 18 : call wrtout([ab_out,std_out],msg,'COLL')
1161 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'xy', flexois(1,iat,2,1,1),flexois(1,iat,2,2,2),flexois(1,iat,2,3,3),&
1162 12 : & flexois(1,iat,2,2,3),flexois(1,iat,2,1,3),flexois(1,iat,2,1,2)
1163 18 : call wrtout([ab_out,std_out],msg,'COLL')
1164 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zy', flexois(3,iat,2,1,1),flexois(3,iat,2,2,2),flexois(3,iat,2,3,3),&
1165 12 : & flexois(3,iat,2,2,3),flexois(3,iat,2,1,3),flexois(3,iat,2,1,2)
1166 18 : call wrtout([ab_out,std_out],msg,'COLL')
1167 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'zx', flexois(3,iat,1,1,1),flexois(3,iat,1,2,2),flexois(3,iat,1,3,3),&
1168 12 : & flexois(3,iat,1,2,3),flexois(3,iat,1,1,3),flexois(3,iat,1,1,2)
1169 18 : call wrtout([ab_out,std_out],msg,'COLL')
1170 6 : write(msg,'(2x,i3,3x,a3,2x,6(f12.6,1x))') iat, 'yx', flexois(2,iat,1,1,1),flexois(2,iat,1,2,2),flexois(2,iat,1,3,3),&
1171 12 : & flexois(2,iat,1,2,3),flexois(2,iat,1,1,3),flexois(2,iat,1,1,2)
1172 20 : call wrtout([ab_out,std_out],msg,'COLL')
1173 : end do
1174 :
1175 2 : write(msg,'(3a)')ch10,' Type-II lattice contribution to flexoelectric tensor (units= nC/m) ',ch10
1176 6 : call wrtout([ab_out,std_out],msg,'COLL')
1177 2 : write(msg,*)' xx yy zz yz xz xy'
1178 6 : call wrtout([ab_out,std_out],msg,'COLL')
1179 14 : do ivar=1,6
1180 12 : elfd=alpha(ivar)
1181 12 : qvecd=beta(ivar)
1182 12 : write(msg,'(3x,a2,6f12.6)') voigt(ivar),lattflexo(elfd,qvecd,1,1),lattflexo(elfd,qvecd,2,2),lattflexo(elfd,qvecd,3,3),&
1183 24 : lattflexo(elfd,qvecd,2,3),lattflexo(elfd,qvecd,1,3),lattflexo(elfd,qvecd,1,2)
1184 38 : call wrtout([ab_out,std_out],msg,'COLL')
1185 : end do
1186 8 : do ivar=4,6
1187 6 : elfd=beta(ivar)
1188 6 : qvecd=alpha(ivar)
1189 6 : write(msg,'(3x,a2,6f12.6)') voigt(ivar+3),lattflexo(elfd,qvecd,1,1),lattflexo(elfd,qvecd,2,2),lattflexo(elfd,qvecd,3,3),&
1190 12 : lattflexo(elfd,qvecd,2,3),lattflexo(elfd,qvecd,1,3),lattflexo(elfd,qvecd,1,2)
1191 20 : call wrtout([ab_out,std_out],msg,'COLL')
1192 : end do
1193 :
1194 : end if
1195 : DBG_EXIT("COLL")
1196 :
1197 2 : end subroutine dtlattflexo
1198 : !!***
1199 :
1200 : !!****f* m_ddb/dm_psinv
1201 : !! NAME
1202 : !! dm_psinv
1203 : !!
1204 : !! FUNCTION
1205 : !! Computes the pseudoinverse of the dynamical matrix (PRB 72, 035105 (2005)).
1206 : !! This piece of code is copied from m_ddb_internalstr.F90.
1207 : !!
1208 : !! INPUTS
1209 : !! asr= if /=0 acustic sume rule is imposed on the dynamical matrix
1210 : !! blkval(2,3,mpert,3,mpert)= 2nd derivatives wrt two atom displacements (at least)
1211 : !! d2asr(2,3,natom,3,natom)=ASR-correction
1212 : !! iout=out file number
1213 : !! mpert=maximum number of ipert
1214 : !! natom=number of atoms in unit cell
1215 : !!
1216 : !! OUTPUT
1217 : !! kmatrix(3*natom,3*natom) = array with the pseudo-inverse of dynamical matrix
1218 : !!
1219 : !! SOURCE
1220 :
1221 2 : subroutine dm_psinv(asr,blkval,d2asr,iout,kmatrix,mpert,natom)
1222 :
1223 : !Arguments -------------------------------
1224 : !scalars
1225 : integer,intent(in) :: asr,iout,mpert,natom
1226 : !arrays
1227 : real(dp),intent(in) :: blkval(2,3,mpert,3,mpert)
1228 : real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
1229 : real(dp),intent(out) :: kmatrix(3*natom,3*natom)
1230 :
1231 : !Local variables -------------------------
1232 : !scalars
1233 : integer :: ii1,ivarA
1234 : integer :: ivarB
1235 : integer :: ier
1236 : character(len=500) :: message
1237 : !arrays
1238 4 : real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
1239 4 : real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
1240 4 : real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
1241 4 : real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom)
1242 4 : real(dp) :: eigval(3*natom-3),eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3)
1243 4 : real(dp) :: eigvecp(2,3*natom,3*natom)
1244 4 : real(dp) :: zhpev1(2,2*3*natom-4)
1245 4 : real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
1246 4 : real(dp) :: d2cart(2,3*natom,3*natom)
1247 :
1248 : ! *********************************************************************
1249 :
1250 : DBG_ENTER("COLL")
1251 :
1252 204 : d2cart(1,:,:) = RESHAPE(blkval(1,1:3,1:natom,1:3,1:natom), SHAPE=[3*natom,3*natom])
1253 200 : d2cart(2,:,:) = zero
1254 :
1255 : !Eventually impose the acoustic sum rule
1256 : !FIXME: this might depend on ifcflag: impose that it is 0 or generalize
1257 2 : call asria_corr(asr,d2asr,d2cart,natom,natom)
1258 : !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
1259 200 : kmatrix(:,:) = d2cart(1,:,:)
1260 :
1261 : !DEBUG
1262 : !write(std_out,'(/,a,/)')'the force constant matrix'
1263 : !do ivarA=1,3*natom
1264 : !write(std_out,'(/)')
1265 : !do ivarB=1,3*natom
1266 : !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA)
1267 : !end do
1268 : !end do
1269 : !ENDDEBUG
1270 :
1271 200 : Nmatr(:,:)=0.0_dp
1272 8 : do ivarB=0,natom-1
1273 28 : do ivarA=0,natom-1
1274 20 : Nmatr(3*ivarA+1,3*ivarB+1) = one
1275 20 : Nmatr(3*ivarA+2,3*ivarB+2) = one
1276 26 : Nmatr(3*ivarA+3,3*ivarB+3) = one
1277 : end do
1278 : end do
1279 :
1280 : !DEBUG
1281 : !do ivarA=1,3*natom
1282 : !write(std_out,'(/)')
1283 : !do ivarB=1,3*natom
1284 : !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA)
1285 : !end do
1286 : !end do
1287 : !ENDDEBUG
1288 :
1289 : !starting the pseudoinverting processes
1290 : !then get the eigenvectors of the big matrix,give values to matrixBp
1291 299 : Bpmatr=0.0_dp
1292 : ii1=1
1293 20 : do ivarA=1,3*natom
1294 119 : do ivarB=1,ivarA
1295 99 : Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
1296 117 : ii1=ii1+1
1297 : end do
1298 : end do
1299 :
1300 : !Bpmatr(2,:) is the imaginary part of the force matrix
1301 : !then call the subroutines CHPEV and ZHPEV to get the eigenvectors
1302 2 : call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
1303 2 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
1304 :
1305 : !DEBUG
1306 : !the eigenval and eigenvec
1307 : !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
1308 : !do ivarA=1,3*natom
1309 : !write(std_out,'(/)')
1310 : !write(std_out,'(es16.6)')eigvalp(ivarA)
1311 : !end do
1312 : !do ivarA=1,3*natom
1313 : !write(std_out,'(/)')
1314 : !do ivarB=1,3*natom
1315 : !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
1316 : !end do
1317 : !end do
1318 : !ENDDEBUG
1319 :
1320 : !Then do the multiplication to get the reduced matrix,in two steps
1321 : !After this the force constant matrix is decouple in two bloks,
1322 : !acoustic and optical ones
1323 6632 : Apmatr(:,:) = MATMUL(TRANSPOSE(eigvecp(1,:,:)), MATMUL(kmatrix, eigvecp(1,:,:)))
1324 :
1325 : !DEBUG
1326 : !the blok diago
1327 : !write(std_out,'(/,a,/)')'matrixAp'
1328 : !do ivarA=1,3*natom
1329 : !write(std_out,'(/)')
1330 : !do ivarB=1,3*natom
1331 : !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
1332 : !end do
1333 : !end do
1334 : !ENDDEBUG
1335 :
1336 : !Check the last three eigenvalues whether too large or not
1337 2 : ivarB=0
1338 8 : do ivarA=3*natom-2,3*natom
1339 8 : if (ABS(Apmatr(ivarA,ivarA))>tol6)then
1340 0 : ivarB=1
1341 : end if
1342 : end do
1343 :
1344 2 : if(ivarB==1)then
1345 0 : write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,&
1346 0 : & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
1347 0 : & ' are too large at Gamma point.',ch10,&
1348 0 : & ' Increase cutoff energy or k-points sampling.',ch10,&
1349 0 : & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
1350 0 : ABI_WARNING(message)
1351 0 : call wrtout(iout,message,'COLL')
1352 : end if
1353 :
1354 : !Give the value of reduced matrix form Apmatr to Amatr
1355 104 : Amatr(:,:) = Apmatr(1:3*natom-3, 1:3*natom-3)
1356 :
1357 : !Now the reduced matrix is in the matrixA, the convert it
1358 : !first give the give the value of matixB from matrixA
1359 : ii1=1
1360 14 : do ivarA=1,3*natom-3
1361 65 : do ivarB=1,ivarA
1362 51 : Bmatr(1,ii1)=Amatr(ivarB,ivarA)
1363 63 : ii1=ii1+1
1364 : end do
1365 : end do
1366 53 : Bmatr(2,:)=0.0_dp
1367 :
1368 : !Call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
1369 2 : call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
1370 2 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
1371 :
1372 : !Check the unstable phonon modes, if the first is negative then print
1373 : !warning message
1374 2 : if(eigval(1)<-1.0*tol8)then
1375 0 : write(message,'(9a)') ch10,&
1376 0 : & ' --- !WARNING',ch10,&
1377 0 : & ' Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,&
1378 0 : & ' The system under calculation is physically unstable.',ch10,&
1379 0 : & ' ---',ch10
1380 0 : call wrtout(std_out,message,'COLL')
1381 : end if
1382 :
1383 : !Do the matrix mutiplication to get pseudoinverse inverse matrix
1384 104 : Amatr(:,:)=0.0_dp
1385 14 : do ivarA=1,3*natom-3
1386 14 : Amatr(ivarA,ivarA)=1.0_dp/eigval(ivarA)
1387 : end do
1388 :
1389 2684 : Cmatr(:,:) = MATMUL(eigvec(1,:,:), MATMUL(Amatr, TRANSPOSE(eigvec(1,:,:))))
1390 :
1391 : !DEBUG
1392 : !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
1393 : !do ivarA=1,3*natom
1394 : !write(std_out,'(/)')
1395 : !do ivarB=1,3*natom
1396 : !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
1397 : !end do
1398 : !end do
1399 : !ENDDEBUG
1400 :
1401 : !So now the inverse of the reduced matrix is in the matrixC
1402 : !now do another mutilplication to get the pseudoinverse of the original
1403 200 : Cpmatr(:,:)=0.0_dp
1404 104 : Cpmatr(1:3*natom-3,1:3*natom-3) = Cmatr(:,:)
1405 :
1406 6632 : Apmatr(:,:) = MATMUL(eigvecp(1,:,:), MATMUL(Cpmatr, TRANSPOSE(eigvecp(1,:,:))))
1407 :
1408 : !Now the inverse is in Cpmatr
1409 200 : kmatrix(:,:)=Apmatr(:,:)
1410 : !transfer the inverse of k-matrix back to the k matrix
1411 : !so now the inverse of k matrix is in the kmatrix
1412 : !ending the part for pseudoinversing the K matrix
1413 :
1414 : DBG_EXIT("COLL")
1415 :
1416 2 : end subroutine dm_psinv
1417 : !!***
1418 :
1419 : !!****f* ABINIT/m_ddb_flexo/ddb_phi1
1420 : !! NAME
1421 : !! ddb_phi1
1422 : !!
1423 : !! FUNCTION
1424 : !! Return the moment of IFCs Phi^(1) and \sum_k Phi^(2)
1425 : !!
1426 : !! INPUTS
1427 : !! ddb<type(ddb_type)>=Long wave 3rd order derivative database.
1428 : !! ddb_version = 8 digit integer giving date. To mantain compatibility with old DDB files.
1429 : !! Crystal<type(crystal_t)>=Crystal structure parameters
1430 : !! filnamddb = name of the ddb file
1431 : !!
1432 : !! OUTPUT
1433 : !! phi1(3,natom,3,natom,3)=First Moment of IFCs Phi^(1)
1434 : !! phi2(3,natom,3,3,3)=Second Moment of IFCs Phi^(2) summed over atoms
1435 : !!
1436 : !! SIDE EFFECTS
1437 : !!
1438 : !! NOTES
1439 : !!
1440 : !! SOURCE
1441 :
1442 1 : subroutine ddb_phi1(ddb,ddb_lw,ddb_version,filnamddb,phi1,phi2,natom)
1443 :
1444 : implicit none
1445 :
1446 : !Arguments ------------------------------------
1447 : !scalars
1448 : integer,intent(in) :: ddb_version, natom
1449 : class(ddb_type),intent(in) :: ddb,ddb_lw
1450 : character(len=fnlen) :: filnamddb
1451 : !arrays
1452 : real(dp), intent(out) :: phi1(3,natom,3,natom,3)
1453 : real(dp), intent(out) :: phi2(3,natom,3,3,3)
1454 :
1455 : !Local variables-------------------------------
1456 : integer,parameter :: cvrsio8=20100401
1457 : integer :: iblok,istrs,strsd1,strsd2,strst,strsd,qvecd
1458 : character(len=500) :: msg
1459 :
1460 : !arrays
1461 : integer :: rfelfd(4),rfphon(4),rfstrs(4),iat,iatd
1462 : integer,parameter :: alpha(6)=(/1,2,3,2,1,1/),beta(6)=(/1,2,3,3,3,2/)
1463 : integer :: rfqvec(4)
1464 : real(dp) :: qphnrm(3),qphon(3,3),fac
1465 2 : real(dp) :: sqrbkt_t1(3,natom,3,3,3)
1466 2 : real(dp) :: d3cart(2,3,ddb%mpert,3,ddb%mpert,3,ddb%mpert)
1467 2 : real(dp) :: phi2tmp(3,ddb%mpert,3,3,3)
1468 :
1469 : ! *********************************************************************
1470 :
1471 : DBG_ENTER("COLL")
1472 :
1473 1 : qphon(:,:)=zero
1474 4 : qphnrm(:)=one
1475 1 : rfphon(:)=0
1476 1 : rfelfd(:)=0
1477 1 : rfstrs(:)=0
1478 1 : rfqvec(:)=0
1479 :
1480 1 : rfphon(1)=1
1481 1 : rfphon(2)=1
1482 1 : rfqvec(3)=1
1483 :
1484 1 : write(msg, '(2a)' ) ch10," Extract the Phi^(1) coeficients from 3DTE"
1485 1 : call wrtout(std_out,msg,'COLL')
1486 1 : call ddb_lw%get_block(iblok,qphon,qphnrm,rfphon,rfelfd,rfstrs,BLKTYP_d3E_lw,rfqvec=rfqvec)
1487 :
1488 1 : if (iblok == 0) then
1489 0 : call wrtout(std_out, " ")
1490 0 : call wrtout(std_out, "--- !WARNING")
1491 0 : call wrtout(std_out, sjoin("- Cannot find Phi^(1) tensor in DDB file:", filnamddb))
1492 0 : call wrtout(std_out, " flexoflag=1 or 3 requires the DDB file to include the corresponding long wave 3rd derivatives")
1493 : end if
1494 :
1495 251339 : d3cart(1,:,:,:,:,:,:) = reshape(ddb_lw%val(1,:,iblok),shape = (/3,ddb%mpert,3,ddb%mpert,3,ddb%mpert/))
1496 251339 : d3cart(2,:,:,:,:,:,:) = reshape(ddb_lw%val(2,:,iblok),shape = (/3,ddb%mpert,3,ddb%mpert,3,ddb%mpert/))
1497 :
1498 : !Define the factors to apply if DDB file has been created with the old version of
1499 : !the longwave driver.
1500 1 : if (ddb_version <= cvrsio8) then
1501 : fac=two
1502 : else
1503 0 : fac=one
1504 : end if
1505 :
1506 : !Extraction of dC/dq and d2C/dqdq
1507 : phi1(:,:,:,:,:) = fac*RESHAPE(d3cart(2,1:3,1:ddb%natom,1:3,1:ddb%natom,1:3,ddb%natom+8), &
1508 634 : & SHAPE=[3,ddb%natom,3,ddb%natom,3])
1509 :
1510 : !Define the factors to apply if DDB file has been created with the old version of
1511 : !the longwave driver. Note that there is a factor -1 to go to the derivatives
1512 1 : if (ddb_version <= cvrsio8) then
1513 : fac=two
1514 : else
1515 0 : fac=-one
1516 : end if
1517 :
1518 472 : phi2 = zero
1519 7 : do istrs=1,6
1520 6 : strsd1=alpha(istrs)
1521 6 : strsd2=beta(istrs)
1522 6 : strst=ddb%natom+3; if (istrs>3) strst=ddb%natom+4
1523 3 : strsd=istrs; if (istrs>3) strsd=istrs-3
1524 25 : do qvecd=1,3
1525 96 : do iat=1,ddb%natom
1526 306 : do iatd=1,3
1527 216 : phi2tmp(iatd,iat,qvecd,strsd1,strsd2)=fac*d3cart(1,iatd,iat,strsd,strst,qvecd,ddb%natom+8)
1528 288 : if (istrs>3) phi2tmp(iatd,iat,qvecd,strsd2,strsd1)=phi2tmp(iatd,iat,qvecd,strsd1,strsd2)
1529 : end do
1530 : end do
1531 : end do
1532 : end do
1533 4 : do qvecd=1,3
1534 13 : do strsd2=1,3
1535 39 : do strsd1=1,3
1536 144 : do iat=1,ddb%natom
1537 : sqrbkt_t1(1:3,iat,strsd1,strsd2,qvecd)=half*(phi2tmp(1:3,iat,qvecd,strsd1,strsd2) + &
1538 459 : & phi2tmp(1:3,iat,strsd2,strsd1,qvecd))
1539 : end do
1540 : end do
1541 : end do
1542 : end do
1543 :
1544 : !Now correct the stress in the square bracketed tesnor tensor
1545 : ! do qvecd=1,3
1546 : ! do strsd2=1,3
1547 : ! do strsd1=1,3
1548 : ! do iatd=1,3
1549 : ! if (iatd==strsd1) then
1550 : ! sqrbkt_t1(iatd,strsd1,strsd2,qvecd)=sqrbkt_t1(iatd,strsd1,strsd2,qvecd) - stress(strsd2,qvecd)
1551 : ! endif
1552 : ! end do
1553 : ! end do
1554 : ! end do
1555 : ! end do
1556 :
1557 : !Now convert back to type-II in order to obtain the frozen ion Lagrange elastic tensor.
1558 472 : phi2(:,:,:,:,:)=zero
1559 4 : do iatd=1,3
1560 13 : do qvecd=1,3
1561 39 : do strsd1=1,3
1562 117 : do strsd2=1,3
1563 432 : do iat=1,ddb%natom
1564 : phi2(iatd,iat,qvecd,strsd1,strsd2)=sqrbkt_t1(iatd,iat,strsd1,qvecd,strsd2) + &
1565 405 : & sqrbkt_t1(iatd,iat,strsd2,strsd1,qvecd)-sqrbkt_t1(iatd,iat,qvecd,strsd2,strsd1)
1566 : end do
1567 : end do
1568 : end do
1569 : end do
1570 : end do
1571 : DBG_EXIT("COLL")
1572 1 : end subroutine ddb_phi1
1573 : !!***
1574 6 : end module m_ddb_flexo
|
