Line data Source code
1 : !!****m* ABINIT/m_ddb_elast
2 : !! NAME
3 : !! m_ddb_elast
4 : !!
5 : !! FUNCTION
6 : !! Elastic properties (clamped-ions and relaxed-ions).
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (XW, DW)
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_elast
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_crystal
28 : use m_ddb
29 : use m_nctk
30 : use netcdf
31 :
32 : use m_fstrings, only : itoa, sjoin
33 : use m_hide_lapack, only : matrginv
34 : use m_dynmat, only : asria_corr
35 : use m_anaddb_dataset, only : anaddb_dataset_type
36 :
37 : implicit none
38 :
39 : private
40 : !!***
41 :
42 : public :: ddb_elast
43 : !!***
44 :
45 : contains
46 : !!***
47 :
48 : !!****f* ABINIT/ddb_elast
49 : !!
50 : !! NAME
51 : !! ddb_elast
52 : !!
53 : !! FUNCTION
54 : !! Get the elastic and compliance tensors, both clamped ion and relaxed ion,
55 : !! under the fixed electric field boundary condition; in which realxed ion
56 : !! tensors can generate two output tensors one is conventional, the other
57 : !! considers the sress correction.
58 : !!
59 : !! INPUTS
60 : !! inp= (derived datatype) contains all the input variables
61 : !! crystal<crystal_t>=Info on crystalline structure.
62 : !! blkval(2,3,mpert,3,mpert,nblok)=
63 : !! second derivatives of total energy with respect to electric fields
64 : !! atom displacements, strain,...... all in cartesian coordinates
65 : !! d2asr= ASR correction to the dynamical matrix at Gamma
66 : !! iblok= bolk number in DDB file
67 : !! iblok_stress= blok number which contain stress tensor
68 : !! instrain=force response internal strain tensor
69 : !! iout=out file number
70 : !! mpert=maximum number of ipert
71 : !! natom=number of atoms in unit cell
72 : !! nblok=number of total bloks in DDB file
73 : !! ncid=NC file handle (open in the caller)
74 : !!
75 : !! OUTPUT
76 : !! elast=relaxed-ion elastic tensor (without stress correction) (6*6) in Voigt notation
77 : !! elast_clamped=clamped-ion elastic tensor (without stress correction) (6*6) in Voigt notation
78 : !! elast_stress=relaxed-ion elastic tensor (with stress correction) (6*6) in Voigt notation
79 : !!
80 : !! NOTES
81 : !! The elastic (compliance) tensors calculated here are under boundary conditions of
82 : !! fixed Electric Field, different from those in ddb_piezo.F90 which are under fixed
83 : !! Displacement Field and incorporate piezoelectric corrections.
84 : !!
85 : !! SOURCE
86 :
87 6 : subroutine ddb_elast(inp,crystal,blkval,compl,compl_clamped,compl_stress,d2asr,&
88 : & elast,elast_clamped,elast_stress,iblok,iblok_stress,&
89 6 : & instrain,iout,mpert,natom,nblok,ncid)
90 :
91 : !Arguments -------------------------------------------
92 : !scalars
93 : integer,intent(in) :: iblok,iblok_stress,iout,mpert,natom,nblok,ncid
94 : !integer,intent(in) :: msize
95 : type(crystal_t),intent(in) :: crystal
96 : type(anaddb_dataset_type),intent(in) :: inp
97 : !arrays
98 : real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),instrain(3*natom,6)
99 : real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
100 : real(dp),intent(out) :: compl(6,6), compl_clamped(6,6),compl_stress(6,6)
101 : real(dp),intent(out) :: elast(6,6), elast_clamped(6,6),elast_stress(6,6)
102 :
103 : !Local variables------------------------------------
104 : !scalars
105 : integer :: ier,ii1,ii2,ipert1,ipert2,ivarA,ivarB,ncerr
106 : real(dp) :: ucvol
107 : logical :: iwrite
108 : character(len=500) :: message
109 : !arrays
110 12 : real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
111 12 : real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
112 12 : real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
113 12 : real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom)
114 12 : real(dp) :: compl_relaxed(6,6),eigval(3*natom-3)
115 12 : real(dp) :: eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3)
116 12 : real(dp) :: eigvecp(2,3*natom,3*natom),elast_relaxed(6,6)
117 12 : real(dp) :: kmatrix(3*natom,3*natom)
118 12 : real(dp) :: new2(6,6),stress(6),zhpev1(2,2*3*natom-4),zhpev1p(2,2*3*natom-1)
119 12 : real(dp) :: zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
120 12 : real(dp) :: d2cart(2,3*natom,3*natom)
121 :
122 : !***************************************************************************
123 6 : compl = zero; compl_clamped = zero; compl_stress = zero
124 6 : elast = zero; elast_clamped = zero; elast_stress = zero
125 :
126 6 : ucvol = crystal%ucvol
127 6 : iwrite = iout > 0
128 :
129 : !extraction of the elastic constants from the blkvals
130 :
131 42 : do ivarA=1,6
132 258 : do ivarB=1,6
133 : ! because the elastic constant is 6*6,
134 : ! so we should judge if the idir is larger than 3 or not
135 216 : if(ivarA>3) then
136 108 : ii1=ivarA-3
137 108 : ipert1=natom+4 !for the shear modulus
138 : else if(ivarA<=3) then
139 108 : ii1=ivarA
140 108 : ipert1=natom+3 !for the diagonal part
141 : end if
142 216 : if(ivarB>3) then
143 108 : ii2=ivarB-3
144 108 : ipert2=natom+4 !for the shear modulus
145 : else if(ivarB<=3) then
146 108 : ii2=ivarB
147 108 : ipert2=natom+3 !for the diagonal part
148 : end if
149 252 : elast(ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok)
150 : end do
151 : end do
152 :
153 : !then consider the volume, because the unit above is in
154 : !Hartree, in fact the elastic constant should be in
155 : !the units of pressure, the energy/volume
156 : !And then transform the unit to si unit using GPa from Hartree/Bohr^3
157 :
158 42 : do ivarA=1,6
159 258 : do ivarB=1,6
160 252 : elast(ivarA,ivarB)=(elast(ivarA,ivarB)/ucvol)*HaBohr3_GPa
161 : end do
162 : end do
163 :
164 : !then should consider the two situations: clamped and relaxed
165 : !ions respectively, give the initial value of elast_clamped
166 6 : elast_clamped(:,:)=elast(:,:)
167 6 : elast_relaxed(:,:)=elast(:,:)
168 :
169 : !then do the matrix mulplication of instrain*K*instrain to get the
170 : !correction of the relaxed ion quantities, in case natom/=1
171 :
172 6 : if( (inp%elaflag==2 .or. inp%elaflag==3 .or. inp%elaflag==4 .or. inp%elaflag==5) .and. natom/=1 )then
173 : ! extracting force matrix at gamma
174 498 : d2cart(1,:,:) = RESHAPE(blkval(1,1:3,1:natom,1:3,1:natom,iblok), (/3*natom,3*natom/))
175 486 : d2cart(2,:,:) = zero
176 :
177 : ! Eventually impose the acoustic sum rule
178 : ! FIXME: this might depend on ifcflag: impose that it is 0 or generalize
179 : !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
180 6 : call asria_corr(inp%asr,d2asr,d2cart,natom,natom)
181 486 : kmatrix = d2cart(1,:,:)
182 :
183 : ! write(std_out,'(/,a,/)')'the k matrix before inverse'
184 : ! do ii1=1,3*natom
185 : ! write(std_out,'(6es16.3)')kmatrix(ii1,1),kmatrix(ii1,2),kmatrix(ii1,3),&
186 : ! & kmatrix(ii1,4),kmatrix(ii1,5),kmatrix(ii1,6)
187 : ! end do
188 :
189 : ! according to formula, invert the kmatrix(3natom,3natom)
190 :
191 : ! NOTE: MJV 13/3/2011 This is just the 3x3 unit matrix copied throughout the dynamical matrix
192 486 : Nmatr(:,:)=zero
193 22 : do ivarB=0,natom-1
194 70 : do ivarA=0,natom-1
195 48 : Nmatr(3*ivarA+1, 3*ivarB+1) = one
196 48 : Nmatr(3*ivarA+2, 3*ivarB+2) = one
197 64 : Nmatr(3*ivarA+3, 3*ivarB+3) = one
198 : end do
199 : end do
200 :
201 : ! The k matrix is not the inverse here - it has not been changed!
202 : ! write(std_out,'(/,a,/)')'the direct inverse of the Kmatrix'
203 : ! do ivarA=1,3*natom
204 : ! write(std_out,'(/)')
205 : ! do ivarB=1,3*natom
206 : ! write(std_out,'(es16.6)')kmatrix(ivarA,ivarB)
207 : ! end do
208 : ! end do
209 :
210 :
211 : ! starting the pseudo-inverse processes then get the eigenvectors of the big matrix,give values to matrixBp
212 : ! Pack the Nmatr matrix in Hermitian form
213 : ii1=1
214 54 : do ivarA=1,3*natom
215 294 : do ivarB=1,ivarA
216 240 : Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
217 288 : ii1=ii1+1
218 : end do
219 : end do
220 246 : Bpmatr(2,:)=zero !the imaginary part of the force matrix
221 : ! then call the subroutines CHPEV and ZHPEV to get the eigenvectors
222 : ! NOTE: MJV there is a huge indeterminacy in this matrix, which has all identical 3x3 block lines
223 : ! this means the orientation of the 0-eigenvalue eigenvectors is kind of random...
224 : ! Is the usage just to get out the translational modes? We know what the eigenvectors look like already!
225 : ! The translational modes are the last 3 with eigenvalue 6
226 : !
227 6 : call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
228 6 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
229 :
230 : ! the eigenval and eigenvec
231 : ! write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
232 : ! do ivarA=1,3*natom
233 : ! write(std_out,'(/)')
234 : ! write(std_out,'(es16.6)')eigvalp(ivarA)
235 : ! end do
236 : ! do ivarA=1,3*natom
237 : ! write(std_out,'(/)')
238 : ! do ivarB=1,3*natom
239 : ! write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
240 : ! end do
241 : ! end do
242 :
243 : ! do the multiplication to get the reduced matrix,in two steps
244 : ! rotate to eigenbasis constructed above to isolate acoustic modes
245 14904 : Apmatr(:,:) = MATMUL(TRANSPOSE(eigvecp(1,:,:)), MATMUL(kmatrix(:,:), eigvecp(1,:,:)))
246 :
247 : ! the blok diagonal parts
248 : ! write(std_out,'(/,a,/)')'Apmatr'
249 : ! do ivarA=1,3*natom
250 : ! write(std_out,'(/)')
251 : ! do ivarB=1,3*natom
252 : ! write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
253 : ! end do
254 : ! end do
255 :
256 : ! check the last three eigenvalues whether too large
257 6 : ivarB=0
258 24 : do ivarA=3*natom-2,3*natom
259 24 : if (ABS(Apmatr(ivarA,ivarA))>tol6) ivarB=1
260 : end do
261 6 : if(ivarB==1)then
262 0 : write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,&
263 0 : & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
264 0 : & ' are too large at Gamma point',ch10,&
265 0 : & ' increase cutoff energy or k-points sampling.',ch10,&
266 0 : & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),&
267 0 : & Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
268 0 : ABI_WARNING(message)
269 0 : call wrtout(iout,message,'COLL')
270 : end if
271 : ! then give the value of reduced matrix form Apmatr to Amatr
272 234 : Amatr(:,:) = Apmatr(1:3*natom-3, 1:3*natom-3)
273 :
274 : ! now the reduced matrix is in the Amatr, the convert it
275 : ! first give the give the value of Bmatr from Amatr
276 : ii1=1
277 36 : do ivarA=1,3*natom-3
278 150 : do ivarB=1,ivarA
279 114 : Bmatr(1,ii1)=Amatr(ivarB,ivarA)
280 144 : ii1=ii1+1
281 : end do
282 : end do
283 120 : Bmatr(2,:)=zero
284 : ! then call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
285 6 : call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
286 6 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
287 :
288 : ! check the unstable phonon modes, if the first is negative then print warning message
289 6 : if(eigval(1)<-1.0*tol8)then
290 0 : write(message,'(a,a,a,a)') ch10,&
291 0 : 'Unstable eigenvalue detected in force constant matrix at Gamma point.',ch10,&
292 0 : 'The system under calculation is physically unstable.'
293 0 : ABI_WARNING(message)
294 0 : call wrtout(iout,message,'COLL')
295 : end if
296 :
297 : ! the do the matrix muplication to get pseudoinverse inverse matrix
298 234 : Cmatr(:,:)=zero
299 36 : do ivarA=1,3*natom-3
300 36 : Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA)
301 : end do
302 5634 : Amatr(:,:) = MATMUL(MATMUL(eigvec(1,:,:), Cmatr(:,:)), TRANSPOSE(eigvec(1,:,:)))
303 :
304 : ! write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
305 : ! do ivarA=1,3*natom
306 : ! write(std_out,'(/)')
307 : ! do ivarB=1,3*natom
308 : ! write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
309 : ! end do
310 : ! end do
311 :
312 : ! so now the inverse of the reduced matrix is in the matrixA
313 : ! now do another mulplication to get the pseudoinverse of the original
314 486 : Cpmatr(:,:)=zero
315 486 : Apmatr(:,:)=zero
316 36 : do ivarA=1,3*natom-3
317 234 : do ivarB=1,3*natom-3
318 228 : Cpmatr(ivarA,ivarB)=Amatr(ivarA,ivarB)
319 : end do
320 : end do
321 :
322 : ! times the eigvecp
323 54 : do ivarA=1,3*natom
324 486 : do ivarB=1,3*natom
325 4800 : do ii1=1,3*natom
326 4752 : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*Cpmatr(ii1,ivarB)
327 : end do
328 : end do
329 : end do
330 486 : Cpmatr(:,:)=zero
331 54 : do ivarA=1,3*natom
332 486 : do ivarB=1,3*natom
333 4800 : do ii1=1,3*natom
334 4752 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1)
335 : end do
336 : end do
337 : end do
338 :
339 : ! now the inverse in in Cpmatr
340 486 : kmatrix(:,:)=Cpmatr(:,:)
341 : ! transfer the inverse of k-matrix back to the k matrix
342 : ! so now the inverse of k matrix is in the kmatrix
343 : ! ending the part for pseudoinversing the K matrix
344 7050 : new2(:,:) = MATMUL(MATMUL(TRANSPOSE(instrain), kmatrix), instrain(:,:))
345 :
346 :
347 : ! then finish the matrix mupl., consider the unit cell volume and the unit change next step
348 258 : new2(:,:)=(new2(:,:)/ucvol)*HaBohr3_GPa
349 :
350 : ! then the relaxed one should be the previous one minus the new2 element
351 258 : elast_relaxed(:,:)=elast_relaxed(:,:)-new2(:,:)
352 : end if
353 : !the above end if end if for elaflag=2 or elafalg=3 or elafalg=4, or elafalg=5 in line 125
354 :
355 : !write(std_out,'(/,a,/)')'debug the unit cell volume'
356 : !write(std_out,'(2es16.6)')ucvol,HaBohr3_GPa
357 :
358 : !then give the initial value of the compl_relaxed(6,6)
359 6 : compl_relaxed(:,:)=elast_relaxed(:,:)
360 :
361 : !*******************************************************************
362 6 : if(inp%elaflag==1.or. inp%elaflag==3)then
363 : ! print out the clamped-ion elastic constants to output file
364 3 : write(message,'(3a)')ch10,' Elastic Tensor (clamped ion) (unit:10^2GP):',ch10
365 3 : call wrtout(std_out,message,'COLL')
366 21 : do ivarA=1,6
367 18 : write(std_out,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,&
368 18 : & elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,&
369 39 : & elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp
370 : end do
371 :
372 3 : call wrtout(iout,message,'COLL')
373 3 : if (iwrite) then
374 21 : do ivarA=1,6
375 18 : write(iout,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,&
376 18 : & elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,&
377 39 : & elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp
378 : end do
379 : end if
380 : end if
381 :
382 6 : if(inp%elaflag==2.or.inp%elaflag==3 .or. inp%elaflag==4.or. inp%elaflag==5)then
383 6 : if(inp%instrflag==0)then
384 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
385 0 : & 'in order to get the elastic tensor(relaxed ion), ',ch10,&
386 0 : & 'one needs information about internal strain ',ch10,&
387 0 : & 'one should set instrflag==1;',ch10,&
388 0 : & 'otherwise the program will continue but give wrong values.'
389 0 : ABI_WARNING(message)
390 0 : call wrtout(iout,message,'COLL')
391 : end if
392 :
393 6 : write(message,'(5a)')ch10,&
394 6 : & ' Elastic Tensor (relaxed ion) (unit:10^2GP):',ch10,&
395 12 : & ' (at fixed electric field boundary condition)',ch10
396 6 : call wrtout(std_out,message,'COLL')
397 42 : do ivarA=1,6
398 36 : write(std_out,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,&
399 36 : & elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,&
400 36 : & elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,&
401 78 : & elast_relaxed(ivarA,6)/100.00_dp
402 : end do
403 :
404 6 : call wrtout(iout,message,'COLL')
405 6 : if (iwrite) then
406 42 : do ivarA=1,6
407 36 : write(iout,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,&
408 36 : & elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,&
409 36 : & elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,&
410 78 : & elast_relaxed(ivarA,6)/100.00_dp
411 : end do
412 : end if
413 : end if
414 :
415 : !then print the corresponding compliances
416 :
417 6 : if(inp%elaflag==1.or.inp%elaflag==3)then
418 : ! compl(:,:)=elast_clamped(:,:) !convert the elastic tensor
419 3 : compl_clamped(:,:)=elast_clamped(:,:)
420 3 : call matrginv(compl_clamped,6,6)
421 3 : write(message,'(a,a,a)')ch10,' Compliance Tensor (clamped ion) (unit: 10^-2GP^-1):',ch10
422 3 : call wrtout(std_out,message,'COLL')
423 :
424 21 : do ivarB=1,6
425 18 : write(std_out,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,&
426 18 : & compl_clamped(ivarB,2)*100.00_dp,&
427 18 : & compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,&
428 18 : & compl_clamped(ivarB,5)*100.00_dp,&
429 39 : & compl_clamped(ivarB,6)*100.00_dp
430 : end do
431 :
432 3 : call wrtout(iout,message,'COLL')
433 :
434 3 : if (iwrite) then
435 21 : do ivarB=1,6
436 18 : write(iout,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,&
437 18 : & compl_clamped(ivarB,2)*100.00_dp,&
438 18 : & compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,&
439 18 : & compl_clamped(ivarB,5)*100.00_dp,&
440 39 : & compl_clamped(ivarB,6)*100.00_dp
441 : end do
442 : end if
443 : end if
444 :
445 6 : if(inp%elaflag==2.or.inp%elaflag==3 .or. inp%elaflag==4 .or. inp%elaflag==5)then
446 : ! compl(:,:)=elast_relaxed(:,:)
447 6 : call matrginv(compl_relaxed,6,6)
448 6 : if(inp%instrflag==0)then
449 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
450 0 : & 'in order to get the compliance tensor(relaxed ion), ',ch10,&
451 0 : & 'one needs information about internal strain ',ch10,&
452 0 : & 'one should set instrflag==1;',ch10,&
453 0 : & 'otherwise the program will continue but give wrong values.'
454 0 : ABI_WARNING(message)
455 0 : call wrtout(iout,message,'COLL')
456 : end if
457 6 : write(message,'(5a)')ch10,&
458 6 : & ' Compliance Tensor (relaxed ion) (unit: 10^-2GP^-1):',ch10,&
459 12 : & ' (at fixed electric field boundary condition)',ch10
460 6 : call wrtout(std_out,message,'COLL')
461 :
462 42 : do ivarB=1,6
463 36 : write(std_out,'(6f12.7)')compl_relaxed(ivarB,1)*100.00_dp,&
464 36 : & compl_relaxed(ivarB,2)*100.00_dp,&
465 36 : & compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,&
466 36 : & compl_relaxed(ivarB,5)*100.00_dp,&
467 78 : & compl_relaxed(ivarB,6)*100.00_dp
468 : end do
469 6 : call wrtout(iout,message,'COLL')
470 :
471 6 : if (iwrite) then
472 42 : do ivarB=1,6
473 36 : write(iout,'(6f12.7)')compl_relaxed(ivarB,1)*100.00,&
474 36 : & compl_relaxed(ivarB,2)*100.00_dp,&
475 36 : & compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,&
476 36 : & compl_relaxed(ivarB,5)*100.00_dp,&
477 78 : & compl_relaxed(ivarB,6)*100.00_dp
478 : end do
479 : end if
480 : end if
481 :
482 : !befor the end , make sure the tensor elast(6,6)
483 : !will have the relaxed ion values and tensor elast_clamped has the clamped ion
484 : !values, and similarily for the corresponding compliance tensors
485 6 : elast_clamped(:,:)=elast(:,:)
486 6 : elast(:,:)=elast_relaxed(:,:)
487 6 : compl(:,:)=compl_relaxed(:,:)
488 :
489 : !begin the part of computing stress corrected elastic tensors
490 6 : if(inp%elaflag==5)then
491 :
492 : ! check the iblok number of first derivative of energy
493 : ! write(std_out,'(/,a,/)')'iblok number at 8:00Pm'
494 : ! write(std_out,'(i)')iblok_stress
495 : ! write(std_out,'(a,f12.7)')'the total energy', blkval(1,1,1)
496 : ! write(std_out,*)'',blkval(1,:,:,:,:,iblok_stress)
497 : ! write(std_out,*)'',blkval(1,:,7,1,1,iblok_stress)
498 :
499 : ! firts give the corect stress values diagonal parts
500 1 : stress(1)=blkval(1,1,natom+3,1,1,iblok_stress)
501 1 : stress(2)=blkval(1,2,natom+3,1,1,iblok_stress)
502 1 : stress(3)=blkval(1,3,natom+3,1,1,iblok_stress)
503 : ! the shear parts
504 1 : stress(4)=blkval(1,1,natom+4,1,1,iblok_stress)
505 1 : stress(5)=blkval(1,2,natom+4,1,1,iblok_stress)
506 1 : stress(6)=blkval(1,3,natom+4,1,1,iblok_stress)
507 : ! then convert the unit from atomic to the GPa unit
508 7 : do ivarA=1,6
509 7 : stress(ivarA)=stress(ivarA)*HaBohr3_GPa
510 : end do
511 : ! give the initial values of elast_stress tensor
512 1 : elast_stress(:,:)=elast_relaxed(:,:)
513 : ! notice that only the first three rows need to be corrected
514 4 : do ivarA=1,3
515 22 : do ivarB=1,6
516 21 : elast_stress(ivarA,ivarB)=elast_stress(ivarA,ivarB)-stress(ivarB)
517 : end do
518 : end do
519 : ! then compute the values of compliance tensor with stress correction
520 1 : compl_stress(:,:)=elast_stress(:,:)
521 1 : call matrginv(compl_stress,6,6)
522 : ! then print out the results of stress corrected elastic and compliance tensors
523 1 : if(inp%instrflag==0)then
524 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
525 0 : & 'In order to get the elastic tensor (relaxed ion with stress correction), ',ch10,&
526 0 : & 'one needs information about internal strain ',ch10,&
527 0 : & 'one should set instrflag==1;',ch10,&
528 0 : & 'otherwise the program will continue but give wrong values.'
529 0 : ABI_WARNING(message)
530 0 : call wrtout(iout,message,'COLL')
531 : end if
532 1 : write(message,'(5a)')ch10,&
533 1 : & ' Elastic Tensor (relaxed ion with stress corrected) (unit:10^2GP)',ch10,&
534 2 : & ' (at fixed electric field boundary condition)',ch10
535 1 : call wrtout(std_out,message,'COLL')
536 1 : call wrtout(iout,message,'COLL')
537 7 : do ivarA=1,6
538 6 : write(std_out,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,&
539 6 : & elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,&
540 13 : & elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp
541 : end do
542 1 : if (iwrite) then
543 7 : do ivarA=1,6
544 6 : write(iout,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,&
545 6 : & elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,&
546 13 : & elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp
547 : end do
548 : end if
549 :
550 : ! then the compliance tensors with stress correction
551 1 : write(message,'(5a)')ch10,&
552 1 : & ' Compliance Tensor (relaxed ion with stress correction) (unit: 10^-2(GP)^-1):',ch10,&
553 2 : & ' (at fixed electric field boundary condition)',ch10
554 1 : call wrtout(std_out,message,'COLL')
555 1 : call wrtout(iout,message,'COLL')
556 7 : do ivarB=1,6
557 6 : write(std_out,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,&
558 6 : & compl_stress(ivarB,2)*100.00_dp,&
559 6 : & compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,&
560 6 : & compl_stress(ivarB,5)*100.00_dp,&
561 13 : & compl_stress(ivarB,6)*100.00_dp
562 : end do
563 1 : if (iwrite) then
564 7 : do ivarB=1,6
565 6 : write(iout,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,&
566 6 : & compl_stress(ivarB,2)*100.00_dp,&
567 6 : & compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,&
568 6 : & compl_stress(ivarB,5)*100.00_dp,&
569 13 : & compl_stress(ivarB,6)*100.00_dp
570 : end do
571 : end if
572 : end if
573 : !end the if 510th line
574 : !end the part of stress corrected elastic and compliance tensors
575 :
576 :
577 : ! Writes the elastic constants tensors (clamped-ion, relaxed-ion with and
578 : ! without stress corrections) to a netCDF file.
579 : !
580 : ! compl=relaxed-ion compliance tensor(without stress correction) (6*6) in Voigt notation
581 : ! compl_clamped=clamped-ion compliance tensor(without stress correction) (6*6) in Voigt notation
582 : ! compl_stress=relaxed-ion compliance tensor(with stress correction) (6*6) in Voigt notation
583 : ! elast=relaxed-ion elastic tensor(without stress correction) (6*6) in Voigt notation
584 : ! elast_clamped=clamped-ion elastic tensor(without stress correction) (6*6) in Voigt notation
585 : ! elast_stress=relaxed-ion elastic tensor(with stress correction) (6*6) in Voigt notation
586 : ! Units are GPa for elastic constants and GPa^-1 for compliance constants
587 :
588 6 : if (ncid /= nctk_noid) then
589 : ! Define dimensions
590 6 : NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
591 :
592 24 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "asr", "elaflag", "instrflag"])
593 6 : NCF_CHECK(ncerr)
594 :
595 : !arrays
596 : ncerr = nctk_def_arrays(ncid, [&
597 : nctkarr_t('internal_strain_tensor', "dp", 'natom3, six'), &
598 : nctkarr_t('compliance_constants_relaxed_ion', "dp", 'six, six'), &
599 : nctkarr_t('compliance_constants_clamped_ion', "dp", 'six, six'), &
600 : nctkarr_t('compliance_constants_relaxed_ion_stress_corrected', "dp", "six, six"), &
601 : nctkarr_t('elastic_constants_relaxed_ion', "dp", 'six, six'), &
602 : nctkarr_t('elastic_constants_clamped_ion', "dp", 'six, six'), &
603 48 : nctkarr_t('elastic_constants_relaxed_ion_stress_corrected', "dp", 'six, six')])
604 6 : NCF_CHECK(ncerr)
605 :
606 : ! Write variables.
607 6 : NCF_CHECK(nctk_set_datamode(ncid))
608 6 : NCF_CHECK(nf90_put_var(ncid, vid('asr'), inp%asr))
609 6 : NCF_CHECK(nf90_put_var(ncid, vid('elaflag'), inp%elaflag))
610 6 : NCF_CHECK(nf90_put_var(ncid, vid('instrflag'), inp%instrflag))
611 6 : NCF_CHECK(nf90_put_var(ncid, vid('internal_strain_tensor'), instrain))
612 6 : NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_relaxed_ion'), compl))
613 6 : NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_clamped_ion'), compl_clamped))
614 6 : NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_relaxed_ion_stress_corrected'), compl_stress))
615 6 : NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_relaxed_ion'), elast))
616 6 : NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_clamped_ion'), elast_clamped))
617 6 : NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_relaxed_ion_stress_corrected'), elast_stress))
618 : end if
619 :
620 : contains
621 :
622 60 : integer function vid(vname)
623 : character(len=*),intent(in) :: vname
624 60 : vid = nctk_idname(ncid, vname)
625 : end function vid
626 :
627 : end subroutine ddb_elast
628 : !!***
629 :
630 18 : end module m_ddb_elast
631 : !!***
|
