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