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