Line data Source code
1 : !!****m* ABINIT/m_ddb_internalstr
2 : !! NAME
3 : !! m_ddb_internalstr
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (XW)
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 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_ddb_internalstr
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_crystal
28 : use m_ddb
29 :
30 : use m_fstrings, only : itoa, sjoin
31 : use m_dynmat, only : asria_corr
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: ddb_internalstr
39 : !!***
40 :
41 : contains
42 : !!***
43 :
44 : !!****f* ABINIT/ddb_internalstr
45 : !!
46 : !! NAME
47 : !! ddb_internalstr
48 : !!
49 : !! FUNCTION
50 : !! Get the insternal strain tensors,both force response and displacement response ones.
51 : !!
52 : !! INPUTS
53 : !! asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
54 : !! blkval(2,3,mpert,3,mpert,nblok)=
55 : !! second derivatives of total energy with respect to electric fields
56 : !! atom displacements,strain,...... all in cartesian coordinates
57 : !! crystal<crystal_t>=Crystalline structure info.
58 : !! iblok= blok number in DDB file
59 : !! iout=out file number
60 : !! mpert=maximum number of ipert
61 : !! msize=Maximum size of dynamical matrices and other perturbations (ddk, dde...)
62 : !! natom=number of atoms in unit cell
63 : !! nblok=number of total bloks in DDB file
64 : !! prt_internalstr=if 2 or higher, print force and displacement internal strain,
65 : !! if 1, print only force internal strain,
66 : !! if 0, do not print internal strain.
67 : !!
68 : !! OUTPUT
69 : !! instrain=force response internal strain tensor
70 : !!
71 : !! NOTES
72 : !! In output of internal strain tensor,column runs from strain1 to
73 : !! strain6(in Voigt notation),row runs from atom1x,atom1y,atom1z,atom2x,.......
74 : !! sum rule is applied on the internal strain tensor
75 : !!
76 : !! SOURCE
77 :
78 35 : subroutine ddb_internalstr(asr,&
79 : !&crystal,&
80 35 : & blkval,&
81 : !&asrq0,&
82 35 : & d2asr,iblok,instrain,iout,mpert,&
83 : !&msize,&
84 : natom,nblok,prt_internalstr)
85 :
86 : !Arguments----------------------------------------------
87 : !scalars
88 : integer,intent(in) :: asr,iblok,iout,mpert,natom,nblok,prt_internalstr
89 : !integer,intent(in) :: msize
90 : !type(crystal_t),intent(in) :: crystal
91 : !type(asrq0_t),intent(inout) :: asrq0
92 : !arrays
93 : real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
94 : real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok)
95 : real(dp),intent(out) :: instrain(3*natom,6)
96 :
97 : !Local variables------------------------------------
98 : !scalars
99 : integer :: idirA,idirB,ier,ii1,ipertA,ipertB,ivarA,ivarB
100 : character(len=500) :: direction,message
101 : !arrays
102 70 : real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
103 70 : real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
104 70 : real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
105 70 : real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom),deviation(3,6)
106 70 : real(dp) :: eigval(3*natom-3),eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3)
107 70 : real(dp) :: eigvecp(2,3*natom,3*natom),instrain_dis(6,3*natom)
108 70 : real(dp) :: kmatrix(3*natom,3*natom),zhpev1(2,2*3*natom-4)
109 70 : real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
110 70 : real(dp) :: d2cart(2,3*natom,3*natom)
111 :
112 : !***************************************************************
113 :
114 : !extract internal strain from DDB matrix data
115 :
116 203 : do ipertA=1,natom
117 707 : do idirA=1,3
118 504 : ivarA=idirA+3*(ipertA-1)
119 3696 : do ivarB=1,6
120 3024 : if(ivarB<=3) then
121 1512 : idirB=ivarB
122 1512 : ipertB=natom+3
123 : ! for the diagonal modulus
124 : else if(ivarB>3) then
125 1512 : idirB=ivarB-3
126 1512 : ipertB=natom+4
127 : ! for the shear modulus
128 : end if
129 3528 : instrain(ivarA,ivarB)=(-1.0_dp)*blkval(1,idirA,ipertA,idirB,ipertB,iblok)
130 : ! write(std_out,'(es16.6)')blkval(1,idirA,ipertA,idirB,ipertB,iblok)
131 : end do
132 : end do
133 : end do
134 : !according to the definition there is a minus sign before the second derivative
135 :
136 : !apply sum rule to the internal strain tensor
137 35 : deviation(:,:)=zero
138 245 : do ivarB=1,6
139 3269 : do ivarA=1,3*natom
140 3024 : if(mod(ivarA,3)==0)then
141 1008 : deviation(1,ivarB)=deviation(1,ivarB)+instrain(ivarA,ivarB)
142 : end if
143 3024 : if(mod(ivarA,3)==1)then
144 1008 : deviation(2,ivarB)=deviation(2,ivarB)+instrain(ivarA,ivarB)
145 : end if
146 3234 : if(mod(ivarA,3)==2)then
147 1008 : deviation(3,ivarB)=deviation(3,ivarB)+instrain(ivarA,ivarB)
148 : end if
149 : end do
150 : end do
151 :
152 245 : do ivarB=1,6
153 3269 : do ivarA=1,3*natom
154 3024 : if(mod(ivarA,3)==0)then
155 1008 : instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(1,ivarB)/natom
156 : end if
157 3024 : if(mod(ivarA,3)==1)then
158 1008 : instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(2,ivarB)/natom
159 : end if
160 3234 : if(mod(ivarA,3)==2)then
161 1008 : instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(3,ivarB)/natom
162 : end if
163 : end do
164 : end do
165 : !ending the sum rule
166 :
167 : !print the force response internal strain constants into the output file
168 35 : if(prt_internalstr>0)then
169 35 : write(message,'(a,a,a,a)')ch10,&
170 70 : & ' Force-response internal strain tensor','(Unit:Hartree/bohr)',ch10
171 35 : call wrtout(std_out,message,'COLL')
172 35 : call wrtout(iout,message,'COLL')
173 :
174 35 : write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
175 70 : & 'strainyy','strainzz','strainyz','strainxz','strainxy'
176 35 : call wrtout(std_out,message,'COLL')
177 539 : do ii1=1,3*natom
178 504 : if(mod(ii1,3)==1)then
179 168 : direction='x'
180 336 : elseif(mod(ii1,3)==2)then
181 168 : direction='y'
182 : elseif(mod(ii1,3)==0)then
183 168 : direction='z'
184 : end if
185 504 : write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,' ',direction,&
186 504 : & instrain(ii1,1),instrain(ii1,2),instrain(ii1,3),&
187 1008 : & instrain(ii1,4),instrain(ii1,5),instrain(ii1,6)
188 539 : call wrtout(std_out,message,'COLL')
189 : end do
190 : endif
191 :
192 : !now write into the ddb output file
193 35 : write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
194 70 : & 'strainyy','strainzz','strainyz','strainxz','strainxy'
195 35 : call wrtout(iout,message,'COLL')
196 539 : do ii1=1,3*natom
197 504 : if(mod(ii1,3)==1)then
198 168 : direction='x'
199 336 : elseif(mod(ii1,3)==2)then
200 168 : direction='y'
201 : elseif(mod(ii1,3)==0)then
202 168 : direction='z'
203 : end if
204 504 : write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,' ',direction,&
205 504 : & instrain(ii1,1),&
206 504 : & instrain(ii1,2),instrain(ii1,3),&
207 1008 : & instrain(ii1,4),instrain(ii1,5),instrain(ii1,6)
208 539 : call wrtout(iout,message,'COLL')
209 : end do
210 :
211 : ! ----------------------------------------------------------------------------------------
212 :
213 : !Try to get the displacement response internal strain tensor
214 : !first need the inverse of force constant matrix
215 10293 : d2cart(1,:,:) = RESHAPE(blkval(1,1:3,1:natom,1:3,1:natom,iblok), SHAPE=[3*natom,3*natom])
216 10223 : d2cart(2,:,:) = zero
217 :
218 : !Eventually impose the acoustic sum rule
219 : !FIXME: this might depend on ifcflag: impose that it is 0 or generalize
220 35 : call asria_corr(asr,d2asr,d2cart,natom,natom)
221 : !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
222 10223 : kmatrix = d2cart(1,:,:)
223 10223 : Apmatr(:,:)=kmatrix(:,:)
224 :
225 : !DEBUG
226 : !write(std_out,'(/,a,/)')'the force constant matrix'
227 : !do ivarA=1,3*natom
228 : !write(std_out,'(/)')
229 : !do ivarB=1,3*natom
230 : !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA)
231 : !end do
232 : !end do
233 : !ENDDEBUG
234 :
235 10223 : Nmatr(:,:)=0.0_dp
236 539 : do ivarA=1,3*natom
237 10223 : do ivarB=1,3*natom
238 9684 : if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0)then
239 1076 : Nmatr(ivarA,ivarB)=one
240 : end if
241 9684 : if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1)then
242 1076 : Nmatr(ivarA,ivarB)=one
243 : end if
244 10188 : if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2)then
245 1076 : Nmatr(ivarA,ivarB)=one
246 : end if
247 : end do
248 : end do
249 :
250 : !DEBUG
251 : !do ivarA=1,3*natom
252 : !write(std_out,'(/)')
253 : !do ivarB=1,3*natom
254 : !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA)
255 : !end do
256 : !end do
257 : !ENDDEBUG
258 :
259 35 : if (natom > 1) then
260 :
261 : !starting the pseudoinverting processes
262 : !then get the eigenvectors of the big matrix,give values to matrixBp
263 15317 : Bpmatr=0.0_dp
264 : ii1=1
265 539 : do ivarA=1,3*natom
266 5633 : do ivarB=1,ivarA
267 5094 : Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
268 5598 : ii1=ii1+1
269 : end do
270 : end do
271 :
272 : !Bpmatr(2,:) is the imaginary part of the force matrix
273 : !then call the subroutines CHPEV and ZHPEV to get the eigenvectors
274 35 : call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
275 35 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
276 :
277 : !DEBUG
278 : !the eigenval and eigenvec
279 : !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
280 : !do ivarA=1,3*natom
281 : !write(std_out,'(/)')
282 : !write(std_out,'(es16.6)')eigvalp(ivarA)
283 : !end do
284 : !do ivarA=1,3*natom
285 : !write(std_out,'(/)')
286 : !do ivarB=1,3*natom
287 : !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
288 : !end do
289 : !end do
290 : !ENDDEBUG
291 :
292 : !Then do the multiplication to get the reduced matrix,in two steps
293 : !After this the force constant matrix is decouple in two bloks,
294 : !acoustic and optical ones
295 10223 : Cpmatr(:,:)=0.0_dp
296 539 : do ivarA=1,3*natom
297 10223 : do ivarB=1,3*natom
298 311832 : do ii1=1,3*natom
299 311328 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*Apmatr(ii1,ivarB)
300 : end do
301 : end do
302 : end do
303 :
304 10223 : Apmatr(:,:)=0.0_dp
305 539 : do ivarA=1,3*natom
306 10223 : do ivarB=1,3*natom
307 311832 : do ii1=1,3*natom
308 311328 : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*eigvecp(1,ii1,ivarB)
309 : end do
310 : end do
311 : end do
312 :
313 : !DEBUG
314 : !the blok diago
315 : !write(std_out,'(/,a,/)')'matrixAp'
316 : !do ivarA=1,3*natom
317 : !write(std_out,'(/)')
318 : !do ivarB=1,3*natom
319 : !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
320 : !end do
321 : !end do
322 : !ENDDEBUG
323 :
324 : !Check the last three eigenvalues whether too large or not
325 : ivarB=0
326 140 : do ivarA=3*natom-2,3*natom
327 140 : if (ABS(Apmatr(ivarA,ivarA))>tol6)then
328 0 : ivarB=1
329 : end if
330 : end do
331 :
332 35 : if(ivarB==1)then
333 0 : write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,&
334 0 : & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
335 0 : & ' are too large at Gamma point.',ch10,&
336 0 : & ' Increase cutoff energy or k-points sampling.',ch10,&
337 0 : & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
338 0 : ABI_WARNING(message)
339 0 : call wrtout(iout,message,'COLL')
340 : end if
341 :
342 : !Give the value of reduced matrix form Apmatr to Amatr
343 434 : do ivarA=1,3*natom-3
344 7409 : do ivarB=1,3*natom-3
345 7374 : Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)
346 : end do
347 : end do
348 :
349 : !Now the reduced matrix is in the matrixA, the convert it
350 : !first give the give the value of matixB from matrixA
351 : ii1=1
352 434 : do ivarA=1,3*natom-3
353 4121 : do ivarB=1,ivarA
354 3687 : Bmatr(1,ii1)=Amatr(ivarB,ivarA)
355 4086 : ii1=ii1+1
356 : end do
357 : end do
358 3722 : Bmatr(2,:)=0.0_dp
359 :
360 : !Call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
361 35 : call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
362 35 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
363 :
364 : !Check the unstable phonon modes, if the first is negative then print
365 : !warning message
366 35 : if(eigval(1)<-1.0*tol8)then
367 8 : write(message,'(9a)') ch10,&
368 8 : & ' --- !WARNING',ch10,&
369 8 : & ' Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,&
370 8 : & ' The system under calculation is physically unstable.',ch10,&
371 16 : & ' ---',ch10
372 8 : call wrtout(std_out,message,'COLL')
373 : end if
374 :
375 : !Do the matrix mutiplication to get pseudoinverse inverse matrix
376 7409 : Cmatr(:,:)=0.0_dp
377 7409 : Amatr(:,:)=0.0_dp
378 434 : do ivarA=1,3*natom-3
379 434 : Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA)
380 : end do
381 :
382 434 : do ivarA=1,3*natom-3
383 7409 : do ivarB=1,3*natom-3
384 234525 : do ii1=1,3*natom-3
385 234126 : Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*Cmatr(ii1,ivarB)
386 : end do
387 : end do
388 : end do
389 :
390 :
391 : !The second multiplication
392 7409 : Cmatr(:,:)=0.0_dp
393 434 : do ivarA=1,3*natom-3
394 7409 : do ivarB=1,3*natom-3
395 234525 : do ii1=1,3*natom-3
396 234126 : Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+ Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1)
397 : end do
398 : end do
399 : end do
400 :
401 : !DEBUG
402 : !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
403 : !do ivarA=1,3*natom
404 : !write(std_out,'(/)')
405 : !do ivarB=1,3*natom
406 : !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
407 : !end do
408 : !end do
409 : !ENDDEBUG
410 :
411 : !So now the inverse of the reduced matrix is in the matrixC
412 : !now do another mutilplication to get the pseudoinverse of the original
413 10223 : Cpmatr(:,:)=0.0_dp
414 10223 : Apmatr(:,:)=0.0_dp
415 434 : do ivarA=1,3*natom-3
416 7409 : do ivarB=1,3*natom-3
417 7374 : Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)
418 : end do
419 : end do
420 :
421 : !Now times the eigvecp
422 539 : do ivarA=1,3*natom
423 10223 : do ivarB=1,3*natom
424 311832 : do ii1=1,3*natom
425 : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*&
426 311328 : & Cpmatr(ii1,ivarB)
427 : end do
428 : end do
429 : end do
430 10223 : Cpmatr(:,:)=0.0_dp
431 539 : do ivarA=1,3*natom
432 10223 : do ivarB=1,3*natom
433 311832 : do ii1=1,3*natom
434 311328 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+ Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1)
435 : end do
436 : end do
437 : end do
438 :
439 : !Now the inverse is in Cpmatr
440 10223 : kmatrix(:,:)=Cpmatr(:,:)
441 : !transfer the inverse of k-matrix back to the k matrix
442 : !so now the inverse of k matrix is in the kmatrix
443 : !ending the part for pseudoinversing the K matrix
444 :
445 : !Now do simple mulplication to obtain the displacement response
446 : !internal strain tensor
447 3563 : instrain_dis(:,:)=0.0_dp
448 245 : do ivarA=1,6
449 3269 : do ivarB=1,3*natom
450 61338 : do ii1=1,3*natom
451 : instrain_dis(ivarA,ivarB)=instrain_dis(ivarA,ivarB)+&
452 61128 : & instrain(ii1,ivarA)*kmatrix(ii1,ivarB)
453 : end do
454 : end do
455 : end do
456 :
457 : else
458 0 : instrain_dis(:,:)=0.0_dp
459 : end if
460 :
461 : !Print out the results
462 35 : if(prt_internalstr>1)then
463 6 : write(message,'(a,a,a,a)')ch10,&
464 12 : & ' Displacement-response internal strain ', 'tensor (Unit:Bohr)',ch10
465 6 : call wrtout(std_out,message,'COLL')
466 6 : call wrtout(iout,message,'COLL')
467 6 : write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
468 12 : & 'strainyy','strainzz','strainyz','strainxz','strainxy'
469 6 : call wrtout(std_out,message,'COLL')
470 6 : call wrtout(iout,message,'COLL')
471 54 : do ivarA=1,3*natom
472 48 : if(mod(ivarA,3)==1)then
473 16 : direction='x'
474 32 : elseif(mod(ivarA,3)==2)then
475 16 : direction='y'
476 : elseif(mod(ivarA,3)==0)then
477 16 : direction='z'
478 : end if
479 48 : write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ivarA-1)/3)+1,' ',direction,&
480 48 : & instrain_dis(1,ivarA),instrain_dis(2,ivarA),&
481 48 : & instrain_dis(3,ivarA),instrain_dis(4,ivarA),instrain_dis(5,ivarA),&
482 96 : & instrain_dis(6,ivarA)
483 48 : call wrtout(std_out,message,'COLL')
484 54 : call wrtout(iout,message,'COLL')
485 : end do
486 : endif
487 :
488 35 : end subroutine ddb_internalstr
489 : !!***
490 :
491 : end module m_ddb_internalstr
492 : !!***
|
