Branch data 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
+ - ]
|
