TABLE OF CONTENTS


ABINIT/ddb_elast [ Functions ]

[ Top ] [ Functions ]

NAME

 ddb_elast

FUNCTION

 Get the elastic and compliance tensors, both clamped ion and relaxed ion,
 under the fixed electric field boundary condition; in which realxed ion
 tensors can generate two output tensors one is conventional, the other
 considers the sress correction.

INPUTS

 inp= (derived datatype) contains all the input variables
 crystal<crystal_t>=Info on crystalline structure.
 blkval(2,3,mpert,3,mpert,nblok)=
   second derivatives of total energy with respect to electric fields
   atom displacements,strain,...... all in cartesian coordinates
 d2asr= ASR correction to the dynamical matrix at Gamma
 iblok= bolk number in DDB file
 iblok_stress= blok number which contain stress tensor
 instrain=force response internal strain tensor
 iout=out file number
 mpert=maximum number of ipert
 natom=number of atoms in unit cell
 nblok=number of total bloks in DDB file
 ncid=NC file handle (open in the caller)

OUTPUT

 elast=relaxed-ion elastic tensor(without stress correction) (6*6) in Voigt notation
 elast_clamped=clamped-ion elastic tensor(without stress correction) (6*6) in Voigt notation
 elast_stress=relaxed-ion elastic tensor(with stress correction) (6*6) in Voigt notation

NOTES

 The elastic (compliance) tensors calculated here are under boundary conditions of
 fixed Electric Field, different from those in ddb_piezo.F90 which are under fixed
 Displacement Field and incorporate piezoelectric corrections.

PARENTS

      anaddb

CHILDREN

      asria_corr,matrginv,wrtout,zhpev

SOURCE

100 subroutine ddb_elast(inp,crystal,blkval,compl,compl_clamped,compl_stress,d2asr,&
101 &            elast,elast_clamped,elast_stress,iblok,iblok_stress,&
102 &            instrain,iout,mpert,natom,nblok,ncid)
103 
104 
105 !This section has been created automatically by the script Abilint (TD).
106 !Do not modify the following lines by hand.
107 #undef ABI_FUNC
108 #define ABI_FUNC 'ddb_elast'
109 !End of the abilint section
110 
111  implicit none
112 
113 !Arguments -------------------------------------------
114 !scalars
115  integer,intent(in) :: iblok,iblok_stress,iout,mpert,natom,nblok,ncid
116 !integer,intent(in) :: msize
117  type(crystal_t),intent(in) :: crystal
118  type(anaddb_dataset_type),intent(in) :: inp
119 !arrays
120  real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),instrain(3*natom,6)
121  real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
122  real(dp),intent(out) :: compl(6,6), compl_clamped(6,6),compl_stress(6,6)
123  real(dp),intent(out) :: elast(6,6), elast_clamped(6,6),elast_stress(6,6)
124 
125 !Local variables------------------------------------
126 !scalars
127  integer :: ier,ii1,ii2,ipert1,ipert2,ivarA,ivarB,ncerr
128  real(dp) :: ucvol
129  logical :: iwrite
130  character(len=500) :: message
131 !arrays
132  real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
133  real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
134  real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
135  real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom)
136  real(dp) :: compl_relaxed(6,6),eigval(3*natom-3)
137  real(dp) :: eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3)
138  real(dp) :: eigvecp(2,3*natom,3*natom),elast_relaxed(6,6)
139  real(dp) :: kmatrix(3*natom,3*natom),new1(6,3*natom)
140  real(dp) :: new2(6,6),stress(6),zhpev1(2,2*3*natom-4),zhpev1p(2,2*3*natom-1)
141  real(dp) :: zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
142  real(dp) :: d2cart(2,3*natom,3*natom)
143 
144 !***************************************************************************
145  compl = zero; compl_clamped = zero; compl_stress = zero
146  elast = zero; elast_clamped = zero; elast_stress = zero
147 
148  ucvol = crystal%ucvol
149  iwrite = iout > 0
150 
151 !extraction of the elastic constants from the blkvals
152 
153  do ivarA=1,6
154    do ivarB=1,6
155 !    because the elastic constant is 6*6,
156 !    so we should judge if the idir is larger than 3 or not
157      if(ivarA>3) then
158        ii1=ivarA-3
159        ipert1=natom+4  !for the shear modulus
160      else if(ivarA<=3) then
161        ii1=ivarA
162        ipert1=natom+3  !for the diagonal part
163      end if
164      if(ivarB>3) then
165        ii2=ivarB-3
166        ipert2=natom+4  !for the shear modulus
167      else if(ivarB<=3) then
168        ii2=ivarB
169        ipert2=natom+3  !for the diagonal part
170      end if
171      elast(ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok)
172    end do
173  end do
174 
175 !then consider the volume, because the unit above is in
176 !Hartree, in fact the elastic constant should be in
177 !the units of presure, the energy/volume
178 !And then transform the unit to si unit using GPa
179 !from Hartree/Bohr^3
180 
181  do ivarA=1,6
182    do ivarB=1,6
183      elast(ivarA,ivarB)=(elast(ivarA,ivarB)/ucvol)*HaBohr3_GPa
184    end do
185  end do
186 
187 !then should consider the two situations:clamped and relaxed
188 !ions respectively,give the initial value of elast_clamped
189  elast_clamped(:,:)=elast(:,:)
190  elast_relaxed(:,:)=elast(:,:)
191 
192 !then do the matrix mulplication of instrain*K*instrain to get the
193 !correction of the relaxed ion quantities, in case natom/=1
194 
195  if( (inp%elaflag==2 .or. inp%elaflag==3&
196 & .or. inp%elaflag==4 .or. inp%elaflag==5) .and. natom/=1 )then
197 !  extracting force matrix at gamma
198    d2cart = zero
199    do ipert1=1,natom
200      do ii1=1,3
201        ivarA=ii1+3*(ipert1-1)
202        do ipert2=1,natom
203          do ii2=1,3
204            ivarB=ii2+3*(ipert2-1)
205            d2cart(1,ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok)
206          end do
207        end do
208      end do
209    end do
210 
211 !  Eventually impose the acoustic sum rule
212 !  FIXME: this might depend on ifcflag: impose that it is 0 or generalize
213    !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
214    call asria_corr(inp%asr,d2asr,d2cart,natom,natom)
215    kmatrix = d2cart(1,:,:)
216 
217 !  DEBUG
218 !  write(std_out,'(/,a,/)')'the k matrix before inverse'
219 !  do ii1=1,3*natom
220 !  write(std_out,'(6es16.3)')kmatrix(ii1,1),kmatrix(ii1,2),kmatrix(ii1,3),&
221 !  & kmatrix(ii1,4),kmatrix(ii1,5),kmatrix(ii1,6)
222 !  end do
223 !  ENDDEBUG
224 
225 !  according to formula, invert the kmatrix(3natom,3natom)
226    Apmatr(:,:)=kmatrix(:,:)
227 
228 !  NOTE: MJV 13/3/2011 This is just the 3x3 unit matrix copied throughout the dynamical matrix
229    Nmatr(:,:)=zero
230    do ivarA=1,3*natom
231      do ivarB=1,3*natom
232        if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0) Nmatr(ivarA,ivarB)=one
233        if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1) Nmatr(ivarA,ivarB)=one
234        if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2) Nmatr(ivarA,ivarB)=one
235      end do
236    end do
237 
238 !  DEBUG
239 !  The k matrix is not the inverse here - it has not been changed!
240 !  write(std_out,'(/,a,/)')'the direct inverse of the Kmatrix'
241 !  do ivarA=1,3*natom
242 !  write(std_out,'(/)')
243 !  do ivarB=1,3*natom
244 !  write(std_out,'(es16.6)')kmatrix(ivarA,ivarB)
245 !  end do
246 !  end do
247 !  ENDDEBUG
248 
249 
250 !  starting the pseudo-inverse processes
251 !  then get the eigenvectors of the big matrix,give values to matrixBp
252 !  Pack the Nmatr matrix in Hermitian form
253    ii1=1
254    do ivarA=1,3*natom
255      do ivarB=1,ivarA
256        Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
257        ii1=ii1+1
258      end do
259    end do
260    Bpmatr(2,:)=zero  !the imaginary part of the force matrix
261 !  then call the subroutines CHPEV and ZHPEV to get the eigenvectors
262 !  NOTE: MJV there is a huge indeterminacy in this matrix, which has all identical 3x3 block lines
263 !  this means the orientation of the 0-eigenvalue eigenvectors is kind of random...
264 !  Is the usage just to get out the translational modes? We know what the eigenvectors look like already!
265 !  The translational modes are the last 3 with eigenvalue 6
266 !
267    call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
268    ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
269 
270 !  DEBUG
271 !  the eigenval and eigenvec
272 !  write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
273 !  do ivarA=1,3*natom
274 !  write(std_out,'(/)')
275 !  write(std_out,'(es16.6)')eigvalp(ivarA)
276 !  end do
277 !  do ivarA=1,3*natom
278 !  write(std_out,'(/)')
279 !  do ivarB=1,3*natom
280 !  write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
281 !  end do
282 !  end do
283 !  ENDDEBUG
284 
285 !  do the multiplication to get the reduced matrix,in two steps
286 !  rotate to eigenbasis constructed above to isolate acoustic modes
287    Cpmatr(:,:)=zero
288    do ivarA=1,3*natom
289      do ivarB=1,3*natom
290        do ii1=1,3*natom
291          Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*&
292 &         Apmatr(ii1,ivarB)
293        end do
294      end do
295    end do
296 
297    Apmatr(:,:)=zero
298    do ivarA=1,3*natom
299      do ivarB=1,3*natom
300        do ii1=1,3*natom
301          Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*&
302 &         eigvecp(1,ii1,ivarB)
303        end do
304      end do
305    end do
306 
307 !  DEBUG
308 !  the blok diagonal parts
309 !  write(std_out,'(/,a,/)')'Apmatr'
310 !  do ivarA=1,3*natom
311 !  write(std_out,'(/)')
312 !  do ivarB=1,3*natom
313 !  write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
314 !  end do
315 !  end do
316 !  ENDDEBUG
317 
318 
319 !  check the last three eigenvalues whether too large
320    ivarB=0
321    do ivarA=3*natom-2,3*natom
322      if (ABS(Apmatr(ivarA,ivarA))>tol6) ivarB=1
323    end do
324    if(ivarB==1)then
325      write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,&
326 &     '  Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
327 &     '  are too large at Gamma point',ch10,&
328 &     '  increase cutoff energy or k-points sampling.',ch10,&
329 &     '  The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),&
330 &     Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
331      MSG_WARNING(message)
332      call wrtout(iout,message,'COLL')
333    end if
334 !  then give the value of reduced matrix form Apmatr to Amatr
335    do ivarA=1,3*natom-3
336      do ivarB=1,3*natom-3
337        Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)
338      end do
339    end do
340 
341 !  now the reduced matrix is in the Amatr, the convert it
342 !  first give the give the value of Bmatr from Amatr
343    ii1=1
344    do ivarA=1,3*natom-3
345      do ivarB=1,ivarA
346        Bmatr(1,ii1)=Amatr(ivarB,ivarA)
347        ii1=ii1+1
348      end do
349    end do
350    Bmatr(2,:)=zero
351 !  then call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
352    call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
353    ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
354 
355 !  check the unstable phonon modes, if the first is negative then print warning message
356    if(eigval(1)<-1.0*tol8)then
357      write(message,'(a,a,a,a)') ch10,&
358 &     'Unstable eigenvalue detected in force constant matrix at Gamma point.',ch10,&
359 &     'The system under calculation is physically unstable.'
360      MSG_WARNING(message)
361      call wrtout(iout,message,'COLL')
362    end if
363 
364 !  the do the matrix muplication to get pseudoinverse inverse matrix
365    Cmatr(:,:)=zero
366    Amatr(:,:)=zero
367    do ivarA=1,3*natom-3
368      Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA)
369    end do
370    do ivarA=1,3*natom-3
371      do ivarB=1,3*natom-3
372        do ii1=1,3*natom-3
373          Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*&
374 &         Cmatr(ii1,ivarB)
375        end do
376      end do
377    end do
378 
379 !  then the second mulplication
380    Cmatr(:,:)=zero
381    do ivarA=1,3*natom-3
382      do ivarB=1,3*natom-3
383        do ii1=1,3*natom-3
384          Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+&
385 &         Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1)
386        end do
387      end do
388    end do
389 
390 !  DEBUG
391 !  write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
392 !  do ivarA=1,3*natom
393 !  write(std_out,'(/)')
394 !  do ivarB=1,3*natom
395 !  write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
396 !  end do
397 !  end do
398 !  ENDDEBUG
399 
400 !  so now the inverse of the reduced matrix is in the matrixC
401 !  now do another mulplication to get the pseudoinverse of the original
402    Cpmatr(:,:)=zero
403    Apmatr(:,:)=zero
404    do ivarA=1,3*natom-3
405      do ivarB=1,3*natom-3
406        Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)
407      end do
408    end do
409 
410 !  times the eigvecp
411    do ivarA=1,3*natom
412      do ivarB=1,3*natom
413        do ii1=1,3*natom
414          Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*&
415 &         Cpmatr(ii1,ivarB)
416        end do
417      end do
418    end do
419    Cpmatr(:,:)=zero
420    do ivarA=1,3*natom
421      do ivarB=1,3*natom
422        do ii1=1,3*natom
423          Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+&
424 &         Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1)
425        end do
426      end do
427    end do
428 
429 !  now the inverse in in Cpmatr
430    kmatrix(:,:)=Cpmatr(:,:)
431 !  transfer the inverse of k-matrix back to the k matrix
432 !  so now the inverse of k matrix is in the kmatrix
433 !  ending the part for pseudoinversing the K matrix
434 !  then do the first matrix mulplication
435    new1(:,:)=zero
436    do ii1=1,6
437      do ii2=1,3*natom
438        do ivarA=1,3*natom
439          new1(ii1,ii2)=new1(ii1,ii2)+instrain(ivarA,ii1)*&
440 &         kmatrix(ivarA,ii2)
441        end do
442      end do
443    end do
444 !  then do the second matrix mulplication, and change the value of kmatrix
445    new2(:,:)=zero
446    do ii1=1,6
447      do ii2=1,6
448        do ivarB=1,3*natom
449          new2(ii1,ii2)=new2(ii1,ii2)+new1(ii1,ivarB)*&
450 &         instrain(ivarB,ii2)
451        end do
452      end do
453    end do
454 !  then finish the matrix mupl., consider the unit cellvolume
455 !  and the unit change next step
456    do ivarA=1,6
457      do ivarB=1,6
458        new2(ivarA,ivarB)=(new2(ivarA,ivarB)/ucvol)*HaBohr3_GPa
459      end do
460    end do
461 !  then the relaxed one should be the previous one minus the new2 element
462    do ivarA=1,6
463      do ivarB=1,6
464        elast_relaxed(ivarA,ivarB)=elast_relaxed(ivarA,ivarB)-&
465 &       new2(ivarA,ivarB)
466      end do
467    end do
468  end if
469 !the above end if end if for elaflag=2 or elafalg=3 or elafalg=4,
470 !or elafalg=5 in line 125
471 
472 !DEBUG
473 !write(std_out,'(/,a,/)')'debug the unit cell volume'
474 !write(std_out,'(2es16.6)')ucvol,HaBohr3_GPa
475 !ENDDEBUG
476 
477 !then give the initial value of the compl_relaxed(6,6)
478  compl_relaxed(:,:)=elast_relaxed(:,:)
479 
480 !*******************************************************************
481  if(inp%elaflag==1.or. inp%elaflag==3)then
482 !  print out the clamped-ion elastic constants to output file
483    write(message,'(3a)')ch10,' Elastic Tensor (clamped ion) (unit:10^2GP):',ch10
484    call wrtout(std_out,message,'COLL')
485    do ivarA=1,6
486      write(std_out,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,&
487 &     elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,&
488 &     elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp
489    end do
490 
491    call wrtout(iout,message,'COLL')
492    if (iwrite) then
493      do ivarA=1,6
494        write(iout,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,&
495 &       elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,&
496 &       elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp
497      end do
498    end if
499  end if
500 
501  if(inp%elaflag==2.or.inp%elaflag==3&
502 & .or. inp%elaflag==4.or. inp%elaflag==5)then
503    if(inp%instrflag==0)then
504      write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
505 &     'in order to get the elastic  tensor(relaxed ion), ',ch10,&
506 &     'one needs information about internal strain ',ch10,&
507 &     'one should set  instrflag==1;',ch10,&
508 &     'otherwise the program will continue but give wrong values.'
509      MSG_WARNING(message)
510      call wrtout(iout,message,'COLL')
511    end if
512 
513    write(message,'(5a)')ch10,&
514 &   ' Elastic Tensor (relaxed ion) (unit:10^2GP):',ch10,&
515 &   '  (at fixed electric field boundary condition)',ch10
516    call wrtout(std_out,message,'COLL')
517    do ivarA=1,6
518      write(std_out,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,&
519 &     elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,&
520 &     elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,&
521 &     elast_relaxed(ivarA,6)/100.00_dp
522    end do
523 
524    call wrtout(iout,message,'COLL')
525    if (iwrite) then
526      do ivarA=1,6
527        write(iout,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,&
528 &       elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,&
529 &       elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,&
530 &       elast_relaxed(ivarA,6)/100.00_dp
531      end do
532    end if
533  end if
534 
535 !then print the corresponding compliances
536 
537  if(inp%elaflag==1.or.inp%elaflag==3)then
538 !  compl(:,:)=elast_clamped(:,:) !convert the elastic tensor
539    compl_clamped(:,:)=elast_clamped(:,:)
540    call matrginv(compl_clamped,6,6)
541    write(message,'(a,a,a)')ch10,' Compliance Tensor (clamped ion) (unit: 10^-2GP^-1):',ch10
542    call wrtout(std_out,message,'COLL')
543 
544    do ivarB=1,6
545      write(std_out,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,&
546 &     compl_clamped(ivarB,2)*100.00_dp,&
547 &     compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,&
548 &     compl_clamped(ivarB,5)*100.00_dp,&
549 &     compl_clamped(ivarB,6)*100.00_dp
550    end do
551 
552    call wrtout(iout,message,'COLL')
553 
554    if (iwrite) then
555      do ivarB=1,6
556        write(iout,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,&
557 &       compl_clamped(ivarB,2)*100.00_dp,&
558 &       compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,&
559 &       compl_clamped(ivarB,5)*100.00_dp,&
560 &       compl_clamped(ivarB,6)*100.00_dp
561      end do
562    end if
563  end if
564 
565  if(inp%elaflag==2.or.inp%elaflag==3&
566 & .or. inp%elaflag==4 .or. inp%elaflag==5)then
567 !  compl(:,:)=elast_relaxed(:,:)
568    call matrginv(compl_relaxed,6,6)
569    if(inp%instrflag==0)then
570      write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
571 &     'in order to get the compliance tensor(relaxed ion), ',ch10,&
572 &     'one needs information about internal strain ',ch10,&
573 &     'one should set  instrflag==1;',ch10,&
574 &     'otherwise the program will continue but give wrong values.'
575      MSG_WARNING(message)
576      call wrtout(iout,message,'COLL')
577    end if
578    write(message,'(5a)')ch10,&
579 &   ' Compliance Tensor (relaxed ion)  (unit: 10^-2GP^-1):',ch10,&
580 &   '  (at fixed electric field boundary condition)',ch10
581    call wrtout(std_out,message,'COLL')
582 
583    do ivarB=1,6
584      write(std_out,'(6f12.7)')compl_relaxed(ivarB,1)*100.00_dp,&
585 &     compl_relaxed(ivarB,2)*100.00_dp,&
586 &     compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,&
587 &     compl_relaxed(ivarB,5)*100.00_dp,&
588 &     compl_relaxed(ivarB,6)*100.00_dp
589    end do
590    call wrtout(iout,message,'COLL')
591 
592    if (iwrite) then
593      do ivarB=1,6
594        write(iout,'(6f12.7)')compl_relaxed(ivarB,1)*100.00,&
595 &       compl_relaxed(ivarB,2)*100.00_dp,&
596 &       compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,&
597 &       compl_relaxed(ivarB,5)*100.00_dp,&
598 &       compl_relaxed(ivarB,6)*100.00_dp
599      end do
600    end if
601  end if
602 
603 !befor the end , make sure the tensor elast(6,6)
604 !will have the relaxed ion values and tensor elast_clamped has the clamped ion
605 !values, and similarily for the corresponding compliance tensors
606  elast_clamped(:,:)=elast(:,:)
607  elast(:,:)=elast_relaxed(:,:)
608  compl(:,:)=compl_relaxed(:,:)
609 
610 !begin the part of computing stress corrected elastic tensors
611  if(inp%elaflag==5)then
612 
613 !  DEBUG
614 !  check the iblok number of first derivative of energy
615 !  write(std_out,'(/,a,/)')'iblok number at 8:00Pm'
616 !  write(std_out,'(i)')iblok_stress
617 !  write(std_out,'(a,f12.7)')'the total energy', blkval(1,1,1)
618 !  write(std_out,*)'',blkval(1,:,:,:,:,iblok_stress)
619 !  write(std_out,*)'',blkval(1,:,7,1,1,iblok_stress)
620 !  ENDDEBUG
621 
622 !  firts give the corect stress values diagonal parts
623    stress(1)=blkval(1,1,natom+3,1,1,iblok_stress)
624    stress(2)=blkval(1,2,natom+3,1,1,iblok_stress)
625    stress(3)=blkval(1,3,natom+3,1,1,iblok_stress)
626 !  the shear parts
627    stress(4)=blkval(1,1,natom+4,1,1,iblok_stress)
628    stress(5)=blkval(1,2,natom+4,1,1,iblok_stress)
629    stress(6)=blkval(1,3,natom+4,1,1,iblok_stress)
630 !  then convert the unit from atomic to the GPa unit
631    do ivarA=1,6
632      stress(ivarA)=stress(ivarA)*HaBohr3_GPa
633    end do
634 !  give the initial values of elast_stress tensor
635    elast_stress(:,:)=elast_relaxed(:,:)
636 !  notice that only the first three rows need to be corrected
637    do ivarA=1,3
638      do ivarB=1,6
639        elast_stress(ivarA,ivarB)=elast_stress(ivarA,ivarB)-stress(ivarB)
640      end do
641    end do
642 !  then compute the values of compliance tensor with stress correction
643    compl_stress(:,:)=elast_stress(:,:)
644    call matrginv(compl_stress,6,6)
645 !  then print out the results of stress corrected elastic and compliance tensors
646    if(inp%instrflag==0)then
647      write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
648 &     'In order to get the elastic tensor (relaxed ion with stress correction), ',ch10,&
649 &     'one needs information about internal strain ',ch10,&
650 &     'one should set  instrflag==1;',ch10,&
651 &     'otherwise the program will continue but give wrong values.'
652      MSG_WARNING(message)
653      call wrtout(iout,message,'COLL')
654    end if
655    write(message,'(5a)')ch10,&
656 &   ' Elastic Tensor (relaxed ion with stress corrected) (unit:10^2GP)',ch10,&
657 &   '  (at fixed electric field boundary condition)',ch10
658    call wrtout(std_out,message,'COLL')
659    call wrtout(iout,message,'COLL')
660    do ivarA=1,6
661      write(std_out,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,&
662 &     elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,&
663 &     elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp
664    end do
665    if (iwrite) then
666      do ivarA=1,6
667        write(iout,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,&
668 &       elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,&
669 &       elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp
670      end do
671    end if
672 
673 !  then the complinace tensors with stress correction
674    write(message,'(5a)')ch10,&
675 &   ' Compliance Tensor (relaxed ion with stress correction) (unit: 10^-2(GP)^-1):',ch10,&
676 &   '  (at fixed electric field boundary condition)',ch10
677    call wrtout(std_out,message,'COLL')
678    call wrtout(iout,message,'COLL')
679    do ivarB=1,6
680      write(std_out,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,&
681 &     compl_stress(ivarB,2)*100.00_dp,&
682 &     compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,&
683 &     compl_stress(ivarB,5)*100.00_dp,&
684 &     compl_stress(ivarB,6)*100.00_dp
685    end do
686    if (iwrite) then
687      do ivarB=1,6
688        write(iout,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,&
689 &       compl_stress(ivarB,2)*100.00_dp,&
690 &       compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,&
691 &       compl_stress(ivarB,5)*100.00_dp,&
692 &       compl_stress(ivarB,6)*100.00_dp
693      end do
694    end if
695  end if
696 !end the if 510th line
697 !end the part of stress corrected elastic and compliance tensors
698 
699 
700  ! Writes the elastic constants tensors (clamped-ion, relaxed-ion with and
701  ! without stress corrections) to a netCDF file.
702  !
703  ! compl=relaxed-ion compliance tensor(without stress correction) (6*6) in Voigt notation
704  ! compl_clamped=clamped-ion compliance tensor(without stress correction) (6*6) in Voigt notation
705  ! compl_stress=relaxed-ion compliance tensor(with stress correction) (6*6) in Voigt notation
706  ! elast=relaxed-ion elastic tensor(without stress correction) (6*6) in Voigt notation
707  ! elast_clamped=clamped-ion elastic tensor(without stress correction) (6*6) in Voigt notation
708  ! elast_stress=relaxed-ion elastic tensor(with stress correction) (6*6) in Voigt notation
709  ! Units are GPa for elastic constants and GPa^-1 for compliance constants
710 
711  if (ncid /= nctk_noid) then
712 #ifdef HAVE_NETCDF
713    ! Define dimensions
714    NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.))
715 
716    ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "asr", "elaflag", "instrflag"])
717    NCF_CHECK(ncerr)
718 
719    !arrays
720    ncerr = nctk_def_arrays(ncid, [&
721      nctkarr_t('internal_strain_tensor', "dp", 'natom3, six'), &
722      nctkarr_t('compliance_constants_relaxed_ion', "dp", 'six, six'), &
723      nctkarr_t('compliance_constants_clamped_ion', "dp", 'six, six'), &
724      nctkarr_t('compliance_constants_relaxed_ion_stress_corrected', "dp", "six, six"), &
725      nctkarr_t('elastic_constants_relaxed_ion', "dp", 'six, six'), &
726      nctkarr_t('elastic_constants_clamped_ion', "dp", 'six, six'), &
727      nctkarr_t('elastic_constants_relaxed_ion_stress_corrected', "dp", 'six, six')])
728    NCF_CHECK(ncerr)
729 
730    ! Write variables.
731    NCF_CHECK(nctk_set_datamode(ncid))
732    NCF_CHECK(nf90_put_var(ncid, vid('asr'), inp%asr))
733    NCF_CHECK(nf90_put_var(ncid, vid('elaflag'), inp%elaflag))
734    NCF_CHECK(nf90_put_var(ncid, vid('instrflag'), inp%instrflag))
735    NCF_CHECK(nf90_put_var(ncid, vid('internal_strain_tensor'), instrain))
736    NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_relaxed_ion'), compl))
737    NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_clamped_ion'), compl_clamped))
738    NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_relaxed_ion_stress_corrected'), compl_stress))
739    NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_relaxed_ion'), elast))
740    NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_clamped_ion'), elast_clamped))
741    NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_relaxed_ion_stress_corrected'), elast_stress))
742 
743 #else
744    MSG_ERROR("Netcdf support not enabled")
745 #endif
746  end if
747 
748 contains
749 
750  integer function vid(vname)
751 
752 
753 !This section has been created automatically by the script Abilint (TD).
754 !Do not modify the following lines by hand.
755 #undef ABI_FUNC
756 #define ABI_FUNC 'vid'
757 !End of the abilint section
758 
759    character(len=*),intent(in) :: vname
760    vid = nctk_idname(ncid, vname)
761  end function vid
762 
763 end subroutine ddb_elast

ABINIT/m_ddb_elast [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddb_elast

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XW, DW)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_ddb_elast
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_crystal
33  use m_ddb
34  use m_nctk
35 #ifdef HAVE_NETCDF
36  use netcdf
37 #endif
38 
39  use m_fstrings,       only : itoa, sjoin
40  use m_hide_lapack,    only : matrginv
41  use m_dynmat,         only : asria_corr
42  use m_anaddb_dataset, only : anaddb_dataset_type
43 
44  implicit none
45 
46  private