Branch data 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 : : !!***
|
