TABLE OF CONTENTS


ABINIT/compute_anharmonics [ Functions ]

[ Top ] [ Functions ]

NAME

 compute_anharmonics

FUNCTION

 Compute strain phonon coupling by finite differences
 Return the effective_potential with the third order

INPUTS

 filenames(17) = path with all name files
 inp <type(multibinit_dtset_type)> = datatype with all the input variables
 comm=MPI communicator

OUTPUT

 eff_pot<type(effective_potential_type)> = effective_potential datatype to be initialized

PARENTS

      multibinit

CHILDREN

      effective_potential_file_read,effective_potential_free
      effective_potential_setelastic3rd,effective_potential_setelastic4th
      effective_potential_setelasticdispcoupling
      effective_potential_setstrainphononcoupling
      effective_potential_writeabiinput,harmonics_terms_applysumrule,ifc_free
      strain_free,strain_get,strain_init,wrtout,xmpi_bcast

SOURCE

 69 subroutine compute_anharmonics(eff_pot,filenames,inp,comm)
 70 
 71  use defs_basis
 72  use m_errors
 73  use m_abicore
 74  use m_xmpi
 75  use m_io_tools, only : open_file
 76 
 77  use m_ifc
 78  use m_anharmonics_terms
 79  use m_effective_potential
 80  use m_effective_potential_file
 81  use m_multibinit_dataset, only : multibinit_dtset_type
 82  use m_strain
 83  use m_fstrings, only : itoa,int2char4,ftoa
 84 
 85 !This section has been created automatically by the script Abilint (TD).
 86 !Do not modify the following lines by hand.
 87 #undef ABI_FUNC
 88 #define ABI_FUNC 'compute_anharmonics'
 89 !End of the abilint section
 90 
 91   implicit none
 92 
 93  !Arguments ------------------------------------
 94  !scalars
 95   integer, intent(in) :: comm
 96   character(len=fnlen),intent(in) :: filenames(17)
 97   type(effective_potential_type),target, intent(inout) :: eff_pot
 98   type(multibinit_dtset_type),intent(in) :: inp
 99  !arrays
100 
101  !Local variables-------------------------------
102  !scalar
103   integer :: ia,ii,ierr,irpt,jj,kk,my_rank,natom
104   integer :: nfile,nrpt,nproc
105   real(dp) :: delta,delta1,delta2
106   character(len=500) :: message
107   character(len=fnlen):: name
108   logical :: files_availables = .True.,has_any_strain = .False.
109   logical :: has_all_strain = .True.
110   logical :: iam_master=.FALSE.
111   integer,parameter :: master=0
112  !arrays
113   integer  :: have_strain(6)
114   real(dp) :: deformation(6,2),elastics3rd(6,6,6)
115   real(dp) :: elastics4th(6,6,6,6),rprimd_def(3,3)
116   type(strain_type) :: strain
117   type(ifc_type) :: phonon_strain(6)
118   logical, allocatable :: file_usable(:)
119   real(dp),allocatable :: elastic_displacement(:,:,:,:)
120   type(effective_potential_type),dimension(:),allocatable :: eff_pots
121   type(strain_type),dimension(:),allocatable :: effpot_strain
122   type(effective_potential_type),pointer :: ref_eff_pot
123 
124  ! *************************************************************************
125 
126   write(message,'(a,(80a),a)') ch10,('=',ii=1,80),ch10
127   call wrtout(ab_out,message,'COLL')
128   call wrtout(std_out,message,'COLL')
129 
130   write(message, '(a,a,a)' )' Compute the third order derivative by finite differences',ch10
131   call wrtout(std_out,message,'COLL')
132   call wrtout(ab_out,message,'COLL')
133 
134   write(message, '(a,a,a)' )' The following files will be used :'
135   call wrtout(std_out,message,'COLL')
136   call wrtout(ab_out,message,'COLL')
137 
138  !==========================================
139  !0)Initialisation of variables:
140 ! Set MPI local varibaless
141   nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
142   iam_master = (my_rank == master)
143 
144  !==========================================
145  !1) Get the list of files
146   nfile = 0
147   jj=6
148   do while (jj < 18)
149     if (filenames(jj)/="") then
150       if(jj==6) nfile = 0
151       write(message, '(a,a)' )'  - ',trim(filenames(jj))
152       call wrtout(std_out,message,'COLL')
153       call wrtout(ab_out,message,'COLL')
154       jj = jj + 1
155       nfile = nfile + 1
156     else
157       exit
158     end if
159   end do
160 
161   if(nfile==0) then
162     write(message,'(a)') '  - No file found -'
163     call wrtout(ab_out,message,'COLL')
164     call wrtout(std_out,message,'COLL')
165   end if
166 
167   write(message,'(a,(80a),a)') ch10,('-',ii=1,80),ch10
168   call wrtout(ab_out,message,'COLL')
169   call wrtout(std_out,message,'COLL')
170 
171  !============================================
172  !2) Read the effectives potential from files"
173  !   - store the reference effective potential
174  !   - Also get the strain
175  !   - perform some checks
176   ABI_DATATYPE_ALLOCATE(eff_pots,(nfile))
177   ABI_DATATYPE_ALLOCATE(effpot_strain,(nfile))
178   ABI_ALLOCATE(file_usable,(nfile))
179 
180   ref_eff_pot => eff_pot
181   file_usable(:) = .True.
182 
183   ii = 1 ! Start at the index 1
184   jj = 6 ! Start at the index 6
185   do while (jj < 18)
186     if (filenames(jj)/="".and.filenames(jj)/="no") then
187       !Read and Intialisation of the effective potential type
188       call effective_potential_file_read(filenames(jj),eff_pots(ii),inp,comm)
189       !Eventualy print the xml file
190 !      if(inp%prt_model==-1.or.inp%prt_model>=3) then
191 !        call int2char4(ii,message)
192 !        name = 'structure_'//trim(itoa(ii-1))//'.xml'
193 !        call isfile(name,'new')
194 !        call effective_potential_writeXML(eff_pots(ii),1,filename=name)
195 !      end if
196 
197       !Fill the eff_pots with the conresponding strain
198       call strain_get(effpot_strain(ii),rprim=eff_pot%crystal%rprimd,&
199 &                     rprim_def=eff_pots(ii)%crystal%rprimd)
200 
201       jj = jj + 1; ii = ii + 1
202 
203       write(message,'(a,(80a))') ch10,('-',ia=1,80)
204       call wrtout(ab_out,message,'COLL')
205       call wrtout(std_out,message,'COLL')
206     else
207       exit
208     end if
209   end do
210 
211   !Do some checks
212   if(iam_master)then
213     do ii=1,size(eff_pots)
214       if (eff_pots(ii)%harmonics_terms%ifcs%nrpt/=ref_eff_pot%harmonics_terms%ifcs%nrpt) then
215         write(message,'(a,I0,a,a,a,a,a,I0,a,a,a,a)' )&
216 &      'the number of cell in reference  (',ref_eff_pot%harmonics_terms%ifcs%nrpt,&
217 &       ') is not equal to the  ',ch10,'the number of cell  in ',trim(filenames(ii+5)),&
218 &      ' (',eff_pots(ii)%harmonics_terms%ifcs%nrpt,')',ch10,'this files cannot be used',ch10
219         MSG_WARNING(message)
220         file_usable(ii) = .False.
221       end if
222       if (eff_pots(ii)%crystal%natom/=ref_eff_pot%crystal%natom) then
223         write(message, '(a,I0,a,a,a,a,a,I0,a,a,a,a)' )&
224 &      'the number of atoms in reference  (',ref_eff_pot%crystal%natom,') is not equal to the  ',ch10,&
225 &      'the number of atoms  in ',trim(filenames(ii+5)),' (',eff_pots(ii)%crystal%natom,')',ch10,&
226 &      'this files cannot be used',ch10
227         MSG_WARNING(message)
228         file_usable(ii) = .False.
229       end if
230       if (eff_pots(ii)%crystal%ntypat/=ref_eff_pot%crystal%ntypat) then
231         write(message, '(a,I0,a,a,a,a,a,I0,a,a,a,a)' )&
232 &      'the number of type of atoms in reference  (',ref_eff_pot%crystal%ntypat,&
233 &       ') is not equal to the  ',&
234 &       ch10,'the number of type of atoms  in ',trim(filenames(ii+5)),&
235 &       ' (',eff_pots(ii)%crystal%ntypat,')',&
236 &       ch10,'this files can not be used',ch10
237         MSG_WARNING(message)
238         file_usable(ii) = .False.
239       end if
240     end do
241   end if
242 
243 ! MPI BROADCAST
244   do ii=1,size(eff_pots)
245     call xmpi_bcast (file_usable(ii), master, comm, ierr)
246   end do
247 
248   if (count((effpot_strain%name=="reference"))>1) then
249     write(message, '(2a)' )&
250 &    ' There is several file corresponding to the reference ',ch10
251     MSG_BUG(message)
252   end if
253 
254   have_strain = 0
255 
256   write(message,'(a)') ' Strains available after reading the files:'
257   call wrtout(ab_out,message,'COLL')
258   call wrtout(std_out,message,'COLL')
259   do ii=1,size(eff_pots)
260     if(effpot_strain(ii)%name /= "".and.file_usable(ii)) then
261       write(message,'(a,a,a,I2,a,(ES10.2),a)')&
262 &       ' A ',trim(effpot_strain(ii)%name),' strain in the direction ',&
263 &       effpot_strain(ii)%direction,' with delta of ',effpot_strain(ii)%delta
264       has_any_strain = .True.
265       call wrtout(ab_out,message,'COLL')
266       call wrtout(std_out,message,'COLL')
267     end if
268   end do
269 
270 
271   if(nfile>1.and.has_any_strain) then
272     write(message,'(a,a,a)') ch10, ' ---analize in more details these files---',ch10
273     call wrtout(ab_out,message,'COLL')
274     call wrtout(std_out,message,'COLL')
275   else
276     write(message,'(a)') '  - No strain found -'
277     call wrtout(ab_out,message,'COLL')
278     call wrtout(std_out,message,'COLL')
279     write(message,'(a,(80a),a)') ch10,('-',ia=1,80),ch10
280     call wrtout(ab_out,message,'COLL')
281     call wrtout(std_out,message,'COLL')
282   end if
283 
284  !First check the strain
285   do ii =1,6
286     jj = 0
287     jj = count(effpot_strain%direction==ii)
288     if(jj>2) then
289       write(message, '(a,I1,a)' )&
290  &    ' There is several file corresponding to strain uniaxial in direction ',ii,ch10
291       MSG_ERROR(message)
292     else
293       name = 'uniaxial'
294       if(ii>=4) name = 'shear'
295       if (jj==1) then
296         write(message, '(a,a,a,I1,a,a)' )&
297 &       ' WARNING: There is only one strain ',trim(name),' in direction ',ii,ch10,&
298 &       '          the finate diferences will not be centering'
299         call wrtout(std_out,message,"COLL")
300         has_all_strain = .False.
301         have_strain(ii)=jj
302       else
303         if(jj==2)then
304           write(message, '(a,a,a,I1,a)' )&
305 &          ' There is two files corresponding to strain ',trim(name),' in direction ',ii,ch10
306           call wrtout(ab_out,message,'COLL')
307           call wrtout(std_out,message,'COLL')
308           have_strain(ii)=jj
309         else
310           write(message, '(a,a,a,I1,a,a)' )&
311 &        ' WARNING: There is no strain ',trim(name),' in direction ',ii,ch10
312           call wrtout(std_out,message,"COLL")
313           has_all_strain = .False.
314           if (inp%strcpling == 2) then
315             do kk = 1,2
316               delta = inp%delta_df
317               if (kk==1) delta = -1 * delta
318               call strain_init(strain,name=name,direction=ii,delta=delta)
319               rprimd_def = matmul(eff_pot%crystal%rprimd,strain%strain)
320               if(kk==1) then
321                 write(message, '(a,a,a,a,a,I1,a,a,a,a)' )&
322 &                 ' if you want to get the correct structure, please run dfpt calculation with',ch10,&
323 &                 ' strain ',trim(name),' in the direction',ii,' with delta=',trim(ftoa(delta)),ch10,&
324 &                 ' The corresponding primitive vectors are:'
325               else
326                 write(message, '(a,a,a,I1,a,a,a,a)' )&
327 &                 ' And a strain ',trim(name),' in the direction',ii,' with delta = ',&
328 &                 trim(ftoa(delta)),ch10,' The corresponding primitive vectors are:'
329               end if
330               call wrtout(ab_out,message,'COLL')
331               call wrtout(std_out,message,'COLL')
332               write(message,'(3(F20.10),a,3(F20.10),a,3(F20.10))')&
333 &               rprimd_def(:,1),ch10, rprimd_def(:,2), ch10,rprimd_def(:,3)
334               call wrtout(ab_out,message,'COLL')
335               call wrtout(std_out,message,'COLL')
336               if(iam_master)then
337                 call effective_potential_writeAbiInput(eff_pot,strain=strain)
338               end if
339               call strain_free(strain)
340             end do
341           end if
342         end if
343       end if
344     end if
345   end do
346 
347 ! check if strain exist
348   if(all(have_strain==0).and.inp%strcpling /= 2) then
349       write(message, '(6a)' )&
350 &    ' WARNING: There is no file corresponding to strain',&
351 &    ' to compute 3rd order derivatives.',ch10,&
352 &    '          In this case the 3rd order derivatives are not set',ch10,&
353 &    '          (add files or set strcpling to 0)'
354       call wrtout(std_out,message,"COLL")
355   end if
356 
357 ! check if the existing strains have opposite deformation
358   deformation = zero
359   do ii=1,6
360     if(have_strain(ii)/=0) then
361       ia = 1
362       do jj=1,size(eff_pots)
363         if (effpot_strain(jj)%direction==ii)then
364           deformation(ii,ia) = effpot_strain(jj)%delta
365           ia = ia + 1
366         end if
367       end do
368       if (have_strain(ii)==2)  then
369         delta1 = deformation(ii,1)
370         delta2 = deformation(ii,2)
371         if (delta1+delta2 > tol15) then
372           write(message, '(a,I1,a,a)' )&
373 &             ' The deformations for strain ',ii,&
374 &             ' are not the opposite',ch10
375           MSG_ERROR(message)
376         end if
377       end if
378     end if
379   end do
380 
381   write(message,'(a,(80a))') ch10,('-',ia=1,80)
382   call wrtout(ab_out,message,'COLL')
383   call wrtout(std_out,message,'COLL')
384 
385   write(message,'(a,a)') ch10, ' After analyzing, the strains available are:'
386   call wrtout(ab_out,message,'COLL')
387   call wrtout(std_out,message,'COLL')
388   if(has_any_strain) then
389     do ii=1,6
390       if(have_strain(ii)/=0) then
391         do jj=1,size(eff_pots)
392           if (effpot_strain(jj)%direction==ii)then
393             write(message,'(a,a,a,I2,a,(ES10.2),a)')&
394 &             ' A ',trim(effpot_strain(jj)%name),' strain in the direction ',&
395 &             effpot_strain(jj)%direction,' with delta of ',effpot_strain(jj)%delta
396             call wrtout(ab_out,message,'COLL')
397             call wrtout(std_out,message,'COLL')
398           end if
399         end do
400       else
401         files_availables = .False.
402       end if
403     end do
404   else
405     files_availables = .False.
406     write(message,'(a)') '  - No strain available -'
407     call wrtout(ab_out,message,'COLL')
408     call wrtout(std_out,message,'COLL')
409   end if
410   write(message,'(a,(80a))') ch10,('-',ia=1,80)
411   call wrtout(ab_out,message,'COLL')
412   call wrtout(std_out,message,'COLL')
413 
414   if(has_all_strain) then
415     write(message,'(3a)') ch10, ' The computation of the third order derivative ',&
416 &    'is possible'
417   else
418     if (inp%strcpling /= 2) then
419       if(ref_eff_pot%has_anharmonicsTerms)then
420         write(message,'(10a)') ch10, ' The computation of the third order derivative ',&
421 &        'is not possible',ch10,' somes files are missing please use strcpling 2 to generate',&
422 &        ' inputs files',ch10,' usable by abinit. The third order derivatives  present in  ',&
423 &        trim(filenames(3)),' will be used'
424       else
425         write(message,'(9a)') ch10, ' The computation of the third order derivative ',&
426 &        'is not possible',ch10,' somes files are missing please use strcpling 2 to generate',&
427 &        ' inputs files',ch10,' usable by abinit. The third order derivative will not be set in',&
428 &        ' the XML file'
429       end if
430     else
431       if(ref_eff_pot%has_anharmonicsTerms)then
432         write(message,'(10a)') ch10, ' The computation of the third order derivative ',&
433 &      'is not possible',ch10,' somes files are missing, the input files usable by abinit have been',&
434 &      ' generate.',ch10,' The third order derivatives present in ',trim(filenames(3)),' will be used'
435       else
436         write(message,'(8a)') ch10, ' The computation of the third order derivative ',&
437 &      'is not possible',ch10,' somes files are missing, the input files usable by abinit have been',&
438 &      ' generate.',ch10,' The third order derivatives will be not set in the XML file'
439       end if
440     end if
441     call wrtout(ab_out,message,'COLL')
442     call wrtout(std_out,message,'COLL')
443   end if
444 
445  !================================================
446  !3) Compute finate differences
447   if(has_all_strain) then
448 
449 !   Allocation of array and set some values
450     nrpt  = ref_eff_pot%harmonics_terms%ifcs%nrpt
451     natom = ref_eff_pot%crystal%natom
452     ABI_ALLOCATE(elastic_displacement,(6,6,3,natom))
453 
454     elastics3rd = zero
455     elastics4th = zero
456 
457     do ii=1,6
458       if(have_strain(ii)/=0) then
459  !      We want the find the index of the perturbation ii in eff_pots(ii)
460  !      And store in delta1 and delta2
461         delta1 = zero
462         delta2 = zero
463         do jj=1,size(eff_pots)
464           if (effpot_strain(jj)%direction==ii.and.(effpot_strain(jj)%direction/=0))then
465             if (abs(delta1)<tol16) then
466               delta1 = jj
467             else
468               delta2 = jj
469             end if
470           end if
471         end do
472         if (abs(delta1)>tol16.and.abs(delta1)>tol16)then
473  !        check if delta1 < delta2, in this case, inverse delta1 and delta2
474           if (effpot_strain(int(delta1))%delta < effpot_strain(int(delta2))%delta) then
475             delta = delta1
476             delta1 = delta2
477             delta2 = delta
478           end if
479 !         Compute strain phonon-coupling
480           phonon_strain(ii)%nrpt =  nrpt
481           ABI_ALLOCATE(phonon_strain(ii)%atmfrc,(3,natom,3,natom,nrpt))
482           ABI_ALLOCATE(phonon_strain(ii)%cell,(3,nrpt))
483           phonon_strain(ii)%atmfrc = zero
484           phonon_strain(ii)%cell =  eff_pots(int(delta1))%harmonics_terms%ifcs%cell
485 
486           do irpt=1,phonon_strain(ii)%nrpt
487             phonon_strain(ii)%atmfrc(:,:,:,:,irpt) =&
488 &           (eff_pots(int(delta1))%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt)&
489 &          - eff_pots(int(delta2))%harmonics_terms%ifcs%atmfrc(:,:,:,:,irpt)) / &
490 &            (2 * abs(effpot_strain(int(delta1))%delta))
491           end do
492 
493           if(inp%asr >= 0) then
494 !           Impose sum rule
495             call harmonics_terms_applySumRule(inp%asr,phonon_strain(ii),&
496 &                                                 eff_pot%crystal%natom)
497           end if
498 
499 !         Compute elastic constants
500           elastics3rd(ii,:,:) = (eff_pots(int(delta1))%harmonics_terms%elastic_constants(:,:)&
501 &          - eff_pots(int(delta2))%harmonics_terms%elastic_constants(:,:)) / &
502 &            (2 * abs(effpot_strain(int(delta1))%delta))
503 
504 !         Compute elastic-displacement coupling
505           elastic_displacement(ii,:,:,:)=(eff_pots(int(delta1))%harmonics_terms%strain_coupling(:,:,:)&
506 &          - eff_pots(int(delta2))%harmonics_terms%strain_coupling(:,:,:)) / &
507 &            (2 * abs(effpot_strain(int(delta1))%delta))
508 
509 !         Compute elastic constants
510           elastics4th(ii,ii,:,:) = (eff_pots(int(delta1))%harmonics_terms%elastic_constants(:,:)&
511 &            - 2*ref_eff_pot%harmonics_terms%elastic_constants(:,:)&
512 &          + eff_pots(int(delta2))%harmonics_terms%elastic_constants(:,:)) / &
513 &            (abs(effpot_strain(int(delta1))%delta)**2)
514         end if
515 
516       end if
517     end do
518 
519 !   Set all the values in the effective potential type
520     call effective_potential_setStrainPhononCoupling(eff_pot,natom,phonon_strain)
521     call effective_potential_setElastic3rd(eff_pot,elastics3rd)
522     call effective_potential_setElastic4th(eff_pot,elastics4th)
523     call effective_potential_setElasticDispCoupling(eff_pot,natom,elastic_displacement)
524 
525 
526 !   Free the phonon-strain coupling array
527     do ii = 1,6
528       call ifc_free(phonon_strain(ii))
529     end do
530     ABI_DEALLOCATE(elastic_displacement)
531 
532     write(message,'(4a)') ch10, ' The computation of the 3rd order elastics constants, ',ch10,&
533 &    ' the phonon-strain coupling and the elastic-displacement coupling is done'
534     call wrtout(ab_out,message,'COLL')
535     call wrtout(std_out,message,'COLL')
536 
537   end if
538 
539 
540  !===============================================
541  !4) Free the array of effective potential
542 
543   do jj=1,nfile
544  !  Free the effective potential type
545     call effective_potential_free(eff_pots(jj))
546   end do
547 
548   ABI_DATATYPE_DEALLOCATE(effpot_strain)
549   ABI_DATATYPE_DEALLOCATE(eff_pots)
550   ABI_DEALLOCATE(file_usable)
551 
552 
553   write(message,'(a,a,a,(80a))') ch10,('=',ii=1,80),ch10
554   call wrtout(ab_out,message,'COLL')
555   call wrtout(std_out,message,'COLL')
556 
557 end subroutine compute_anharmonics

ABINIT/m_compute_anharmonics [ Modules ]

[ Top ] [ Modules ]

NAME

  m_compute_anharmonics

FUNCTION

COPYRIGHT

  Copyright (C) 2008-2018 ABINIT group ()
  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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_compute_anharmonics
27 
28  implicit none
29 
30  private