TABLE OF CONTENTS


ABINIT/ddb_internalstr [ Functions ]

[ Top ] [ Functions ]

NAME

 ddb_internalstr

FUNCTION

 Get the insternal strain tensors,both force response and displacement response ones.

INPUTS

 asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file.
 blkval(2,3,mpert,3,mpert,nblok)=
   second derivatives of total energy with respect to electric fields
   atom displacements,strain,...... all in cartesian coordinates
 crystal<crystal_t>=Crystalline structure info.
 iblok= blok number in DDB file
 iout=out file number
 mpert=maximum number of ipert
 msize=Maximum size of dynamical matrices and other perturbations (ddk, dde...)
 natom=number of atoms in unit cell
 nblok=number of total bloks in DDB file
 prt_internalstr=if 2 or higher, print force and displacement internal strain, 
                 if 1, print only force internal strain, 
                 if 0, do not print internal strain. 

OUTPUT

 instrain=force response internal strain tensor

NOTES

 In output of internal strain tensor,column runs from strain1 to
 strain6(in Voigt notation),row runs from atom1x,atom1y,atom1z,atom2x,.......
 sum rule is applied on the internal strain tensor

PARENTS

      anaddb,m_effective_potential_file

CHILDREN

      asria_corr,wrtout,zhpev

SOURCE

 89 subroutine ddb_internalstr(asr,&
 90 !&crystal,&
 91 & blkval,&
 92 !&asrq0,&
 93 & d2asr,iblok,instrain,iout,mpert,&
 94 !&msize,&
 95 natom,nblok,prt_internalstr)
 96 
 97 
 98 !This section has been created automatically by the script Abilint (TD).
 99 !Do not modify the following lines by hand.
100 #undef ABI_FUNC
101 #define ABI_FUNC 'ddb_internalstr'
102 !End of the abilint section
103 
104  implicit none
105 
106 !Arguments----------------------------------------------
107 !scalars
108  integer,intent(in) :: asr,iblok,iout,mpert,natom,nblok,prt_internalstr
109 !integer,intent(in) :: msize
110 !type(crystal_t),intent(in) :: crystal
111 !type(asrq0_t),intent(inout) :: asrq0
112 !arrays
113  real(dp),intent(in) :: d2asr(2,3,natom,3,natom)
114  real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok)
115  real(dp),intent(out) :: instrain(3*natom,6)
116 
117 !Local variables------------------------------------
118 !scalars
119  integer :: idirA,idirB,ier,ii1,ii2,ipertA,ipertB,ivarA,ivarB
120  character(len=500) :: direction,message
121 !arrays
122  real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
123  real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
124  real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
125  real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom),deviation(3,6)
126  real(dp) :: eigval(3*natom-3),eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3)
127  real(dp) :: eigvecp(2,3*natom,3*natom),instrain_dis(6,3*natom)
128  real(dp) :: kmatrix(3*natom,3*natom),zhpev1(2,2*3*natom-4)
129  real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
130  real(dp) :: d2cart(2,3*natom,3*natom)
131 
132 !***************************************************************
133 
134 !extract internal strain from DDB matrix data
135 
136  do ipertA=1,natom
137    do idirA=1,3
138      ivarA=idirA+3*(ipertA-1)
139      do ivarB=1,6
140        if(ivarB<=3) then
141          idirB=ivarB
142          ipertB=natom+3
143 !        for the diagonal modulus
144        else if(ivarB>3) then
145          idirB=ivarB-3
146          ipertB=natom+4
147 !        for the shear modulus
148        end if
149        instrain(ivarA,ivarB)=(-1.0_dp)*blkval(1,idirA,ipertA,idirB,ipertB,iblok)
150 !      write(std_out,'(es16.6)')blkval(1,idirA,ipertA,idirB,ipertB,iblok)
151      end do
152    end do
153  end do
154 !according to the definition there is a minus sign before the second derivative
155 
156 !apply sum rule to the internal strain tensor
157  deviation(:,:)=zero
158  do ivarB=1,6
159    do ivarA=1,3*natom
160      if(mod(ivarA,3)==0)then
161        deviation(1,ivarB)=deviation(1,ivarB)+instrain(ivarA,ivarB)
162      end if
163      if(mod(ivarA,3)==1)then
164        deviation(2,ivarB)=deviation(2,ivarB)+instrain(ivarA,ivarB)
165      end if
166      if(mod(ivarA,3)==2)then
167        deviation(3,ivarB)=deviation(3,ivarB)+instrain(ivarA,ivarB)
168      end if
169    end do
170  end do
171 
172  do ivarB=1,6
173    do ivarA=1,3*natom
174      if(mod(ivarA,3)==0)then
175        instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(1,ivarB)/natom
176      end if
177      if(mod(ivarA,3)==1)then
178        instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(2,ivarB)/natom
179      end if
180      if(mod(ivarA,3)==2)then
181        instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(3,ivarB)/natom
182      end if
183    end do
184  end do
185 !ending the sum rule
186 
187 !print the force response internal strain constants into the output file
188  if(prt_internalstr>0)then
189    write(message,'(a,a,a,a)')ch10,&
190 &   ' Force-response internal strain tensor','(Unit:Hartree/bohr)',ch10
191    call wrtout(std_out,message,'COLL')
192    call wrtout(iout,message,'COLL')
193 
194    write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
195 &   'strainyy','strainzz','strainyz','strainxz','strainxy'
196    call wrtout(std_out,message,'COLL')
197    do ii1=1,3*natom
198      if(mod(ii1,3)==1)then
199        direction='x'
200      elseif(mod(ii1,3)==2)then
201        direction='y'
202      elseif(mod(ii1,3)==0)then
203        direction='z'
204      end if
205      write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,'  ',direction,&
206 &     instrain(ii1,1),instrain(ii1,2),instrain(ii1,3),&
207 &     instrain(ii1,4),instrain(ii1,5),instrain(ii1,6)
208      call wrtout(std_out,message,'COLL')
209    end do
210  endif
211 
212 !now write into the ddb output file
213  write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
214 & 'strainyy','strainzz','strainyz','strainxz','strainxy'
215  call wrtout(iout,message,'COLL')
216  do ii1=1,3*natom
217    if(mod(ii1,3)==1)then
218      direction='x'
219    elseif(mod(ii1,3)==2)then
220      direction='y'
221    elseif(mod(ii1,3)==0)then
222      direction='z'
223    end if
224    write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,'  ',direction,&
225 &   instrain(ii1,1),&
226 &   instrain(ii1,2),instrain(ii1,3),&
227 &   instrain(ii1,4),instrain(ii1,5),instrain(ii1,6)
228    call wrtout(iout,message,'COLL')
229  end do
230 
231 ! ----------------------------------------------------------------------------------------
232 
233 !Try to get the displacement response internal strain tensor
234 !first need the inverse of force constant matrix
235  d2cart = zero
236  do ipertA=1,natom
237    do ii1=1,3
238      ivarA=ii1+3*(ipertA-1)
239      do ipertB=1,natom
240        do ii2=1,3
241          ivarB=ii2+3*(ipertB-1)
242          d2cart(1,ivarA,ivarB)=blkval(1,ii1,ipertA,ii2,ipertB,iblok)
243        end do
244      end do
245    end do
246  end do
247 
248 !Eventually impose the acoustic sum rule
249 !FIXME: this might depend on ifcflag: impose that it is 0 or generalize
250  call asria_corr(asr,d2asr,d2cart,natom,natom)
251  !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart)
252  kmatrix = d2cart(1,:,:)
253  Apmatr(:,:)=kmatrix(:,:)
254 
255 !DEBUG
256 !write(std_out,'(/,a,/)')'the force constant matrix'
257 !do ivarA=1,3*natom
258 !write(std_out,'(/)')
259 !do ivarB=1,3*natom
260 !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA)
261 !end do
262 !end do
263 !ENDDEBUG
264 
265  Nmatr(:,:)=0.0_dp
266  do ivarA=1,3*natom
267    do ivarB=1,3*natom
268      if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0)then
269        Nmatr(ivarA,ivarB)=one
270      end if
271      if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1)then
272        Nmatr(ivarA,ivarB)=one
273      end if
274      if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2)then
275        Nmatr(ivarA,ivarB)=one
276      end if
277    end do
278  end do
279 
280 !DEBUG
281 !do ivarA=1,3*natom
282 !write(std_out,'(/)')
283 !do ivarB=1,3*natom
284 !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA)
285 !end do
286 !end do
287 !ENDDEBUG
288 
289 !starting the pseudoinverting processes
290 !then get the eigenvectors of the big matrix,give values to matrixBp
291  Bpmatr=0.0_dp
292  ii1=1
293  do ivarA=1,3*natom
294    do ivarB=1,ivarA
295      Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
296      ii1=ii1+1
297    end do
298  end do
299 
300 !Bpmatr(2,:) is the imaginary part of the force matrix
301 !then call the subroutines CHPEV and ZHPEV to get the eigenvectors
302  call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
303  ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
304 
305 !DEBUG
306 !the eigenval and eigenvec
307 !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
308 !do ivarA=1,3*natom
309 !write(std_out,'(/)')
310 !write(std_out,'(es16.6)')eigvalp(ivarA)
311 !end do
312 !do ivarA=1,3*natom
313 !write(std_out,'(/)')
314 !do ivarB=1,3*natom
315 !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
316 !end do
317 !end do
318 !ENDDEBUG
319 
320 !Then do the multiplication to get the reduced matrix,in two steps
321 !After this the force constant matrix is decouple in two bloks,
322 !acoustic and optical ones
323  Cpmatr(:,:)=0.0_dp
324  do ivarA=1,3*natom
325    do ivarB=1,3*natom
326      do ii1=1,3*natom
327        Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*Apmatr(ii1,ivarB)
328      end do
329    end do
330  end do
331 
332  Apmatr(:,:)=0.0_dp
333  do ivarA=1,3*natom
334    do ivarB=1,3*natom
335      do ii1=1,3*natom
336        Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*eigvecp(1,ii1,ivarB)
337      end do
338    end do
339  end do
340 
341 !DEBUG
342 !the blok diago
343 !write(std_out,'(/,a,/)')'matrixAp'
344 !do ivarA=1,3*natom
345 !write(std_out,'(/)')
346 !do ivarB=1,3*natom
347 !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
348 !end do
349 !end do
350 !ENDDEBUG
351 
352 !Check the last three eigenvalues whether too large or not
353  ivarB=0
354  do ivarA=3*natom-2,3*natom
355    if (ABS(Apmatr(ivarA,ivarA))>tol6)then
356      ivarB=1
357    end if
358  end do
359 
360  if(ivarB==1)then
361    write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,&
362 &   '  Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
363 &   '  are too large at Gamma point.',ch10,&
364 &   '  Increase cutoff energy or k-points sampling.',ch10,&
365 &   '  The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
366    MSG_WARNING(message)
367    call wrtout(iout,message,'COLL')
368  end if
369 
370 !Give the value of reduced matrix form Apmatr to Amatr
371  do ivarA=1,3*natom-3
372    do ivarB=1,3*natom-3
373      Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)
374    end do
375  end do
376 
377 !Now the reduced matrix is in the matrixA, the convert it
378 !first give the give the value of matixB from matrixA
379  ii1=1
380  do ivarA=1,3*natom-3
381    do ivarB=1,ivarA
382      Bmatr(1,ii1)=Amatr(ivarB,ivarA)
383      ii1=ii1+1
384    end do
385  end do
386  Bmatr(2,:)=0.0_dp
387 
388 !Call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
389  call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
390  ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
391 
392 !Check the unstable phonon modes, if the first is negative then print
393 !warning message
394  if(eigval(1)<-1.0*tol8)then
395    write(message,'(9a)') ch10,&
396 &   ' --- !WARNING',ch10,&
397 &   '     Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,&
398 &   '     The system under calculation is physically unstable.',ch10,&
399 &   ' ---',ch10
400    call wrtout(std_out,message,'COLL')
401  end if
402 
403 !Do the matrix mutiplication to get pseudoinverse inverse matrix
404  Cmatr(:,:)=0.0_dp
405  Amatr(:,:)=0.0_dp
406  do ivarA=1,3*natom-3
407    Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA)
408  end do
409 
410  do ivarA=1,3*natom-3
411    do ivarB=1,3*natom-3
412      do ii1=1,3*natom-3
413        Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*Cmatr(ii1,ivarB)
414      end do
415    end do
416  end do
417 
418 
419 !The second multiplication
420  Cmatr(:,:)=0.0_dp
421  do ivarA=1,3*natom-3
422    do ivarB=1,3*natom-3
423      do ii1=1,3*natom-3
424        Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+ Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1)
425      end do
426    end do
427  end do
428 
429 !DEBUG
430 !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
431 !do ivarA=1,3*natom
432 !write(std_out,'(/)')
433 !do ivarB=1,3*natom
434 !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
435 !end do
436 !end do
437 !ENDDEBUG
438 
439 !So now the inverse of the reduced matrix is in the matrixC
440 !now do another mutilplication to get the pseudoinverse of the original
441  Cpmatr(:,:)=0.0_dp
442  Apmatr(:,:)=0.0_dp
443  do ivarA=1,3*natom-3
444    do ivarB=1,3*natom-3
445      Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)
446    end do
447  end do
448 
449 !Now times the eigvecp
450  do ivarA=1,3*natom
451    do ivarB=1,3*natom
452      do ii1=1,3*natom
453        Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*&
454 &       Cpmatr(ii1,ivarB)
455      end do
456    end do
457  end do
458  Cpmatr(:,:)=0.0_dp
459  do ivarA=1,3*natom
460    do ivarB=1,3*natom
461      do ii1=1,3*natom
462        Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+ Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1)
463      end do
464    end do
465  end do
466 
467 !Now the inverse is in Cpmatr
468  kmatrix(:,:)=Cpmatr(:,:)
469 !transfer the inverse of k-matrix back to the k matrix
470 !so now the inverse of k matrix is in the kmatrix
471 !ending the part for pseudoinversing the K matrix
472 
473 !Now do simple mulplication to obtain the displacement response
474 !internal strain tensor
475 
476  instrain_dis(:,:)=0.0_dp
477  do ivarA=1,6
478    do ivarB=1,3*natom
479      do ii1=1,3*natom
480        instrain_dis(ivarA,ivarB)=instrain_dis(ivarA,ivarB)+&
481 &       instrain(ii1,ivarA)*kmatrix(ii1,ivarB)
482      end do
483    end do
484  end do
485 
486 !Print out the results
487  if(prt_internalstr>1)then
488    write(message,'(a,a,a,a)')ch10,&
489 &   ' Displacement-response internal strain ', 'tensor (Unit:Bohr)',ch10
490    call wrtout(std_out,message,'COLL')
491    call wrtout(iout,message,'COLL')
492    write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',&
493 &   'strainyy','strainzz','strainyz','strainxz','strainxy'
494    call wrtout(std_out,message,'COLL')
495    call wrtout(iout,message,'COLL')
496    do ivarA=1,3*natom
497      if(mod(ivarA,3)==1)then
498        direction='x'
499      elseif(mod(ivarA,3)==2)then
500        direction='y'
501      elseif(mod(ivarA,3)==0)then
502        direction='z'
503      end if
504      write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ivarA-1)/3)+1,'  ',direction,&
505 &     instrain_dis(1,ivarA),instrain_dis(2,ivarA),&
506 &     instrain_dis(3,ivarA),instrain_dis(4,ivarA),instrain_dis(5,ivarA),&
507 &     instrain_dis(6,ivarA)
508      call wrtout(std_out,message,'COLL')
509      call wrtout(iout,message,'COLL')
510    end do
511  endif
512 
513 end subroutine ddb_internalstr

ABINIT/m_ddb_internalstr [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddb_internalstr

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XW)
  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_internalstr
28 
29  use defs_basis
30  use m_abicore
31  use m_errors
32  use m_crystal
33  use m_ddb
34 
35  use m_fstrings,     only : itoa, sjoin
36  use m_dynmat,       only : asria_corr
37 
38  implicit none
39 
40  private