Branch data Line data Source code
1 : : !!****m* ABINIT/m_ddb_diel
2 : : !! NAME
3 : : !! m_ddb_diel
4 : : !!
5 : : !! FUNCTION
6 : : !! This module provides routines for the calculation of the dielectric constant (anaddb)
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1999-2026 ABINIT group (XG,XW, MVeithen, EB)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_ddb_diel
23 : :
24 : : use defs_basis
25 : : use m_errors
26 : : use m_xmpi
27 : : use m_abicore
28 : : use m_ddb
29 : : use m_nctk
30 : : use netcdf
31 : :
32 : : use m_anaddb_dataset, only : anaddb_dataset_type
33 : : use m_crystal, only : crystal_t
34 : :
35 : : implicit none
36 : :
37 : : private
38 : : !!***
39 : :
40 : : public :: ddb_diel
41 : : public :: alignph
42 : : !!***
43 : :
44 : : contains
45 : : !!***
46 : :
47 : : !!****f* ABINIT/ddb_diel
48 : : !!
49 : : !! NAME
50 : : !! ddb_diel
51 : : !!
52 : : !! FUNCTION
53 : : !! Get the frequency-dependent dielectric matrix, as well as the
54 : : !! oscillator strengths and mode effective charges,
55 : : !! and reflectivities (without damping)
56 : : !! See the definitions Eq.(53-54) in PRB55, 10355 (1997) [[cite:Gonze1997a]].
57 : : !!
58 : : !! INPUTS
59 : : !! amu(ntypat)=mass of the atoms (atomic mass unit)
60 : : !! anaddb_dtset= (derived datatype) contains all the input variables
61 : : !! matrix (diagonal in the atoms)
62 : : !! displ(2,3*natom,3*natom)=
63 : : !! the displacements of atoms in cartesian coordinates.
64 : : !! The first index means either the real or the imaginary part,
65 : : !! The second index runs on the direction and the atoms displaced
66 : : !! The third index runs on the modes.
67 : : !! d2cart(2,3,mpert,3,mpert)=dynamical matrix, effective charges, dielectric tensor,... all in cartesian coordinates
68 : : !! iout=unit number for outputs
69 : : !! lst(3*nph2l)=log. of product of frequencies**2, needed to calculate
70 : : !! the generalized Lyddane-Sachs-Teller relation at zero frequency
71 : : !! mpert =maximum number of ipert
72 : : !! natom=number of atoms in unit cell
73 : : !! nph2l=input variable from anaddb_dtset, needed to dimension lst
74 : : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
75 : : !! matrix eigenvalues, except if these are negative, and in this
76 : : !! case, give minus the square root of the absolute value
77 : : !! of the matrix eigenvalues). Hartree units.
78 : : !! comm=MPI communicator.
79 : : !! ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.
80 : : !!
81 : : !! OUTPUT
82 : : !! fact_oscstr(2,3,3*natom)=oscillator strengths for the different eigenmodes,
83 : : !! for different direction of the electric field;
84 : : !! dielt_rlx(3,3) relaxed ion (zero frequency) dielectric tensor.
85 : : !!
86 : : !! NOTES
87 : : !! 1. The phonon frequencies phfrq should correspond to the
88 : : !! wavevector at Gamma, without any non-analyticities.
89 : : !! 2. Should clean for no imaginary part ...
90 : : !! This routine should be used only by one processor.
91 : : !! 3. frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
92 : : !! mode effective charges for the different eigenmodes,
93 : : !! for different direction of the electric field
94 : : !!
95 : : !! SOURCE
96 : :
97 : 26 : subroutine ddb_diel(Crystal,amu,anaddb_dtset,dielt_rlx,displ,d2cart,epsinf,fact_oscstr,&
98 : 26 : & iout,lst,mpert,natom,nph2l,phfrq,comm,ncid)
99 : :
100 : : !Arguments -------------------------------
101 : : !scalars
102 : : integer,intent(in) :: iout,mpert,natom,nph2l,comm,ncid
103 : : type(crystal_t),intent(in) :: Crystal
104 : : type(anaddb_dataset_type),intent(in) :: anaddb_dtset
105 : : real(dp),intent(in) :: lst(nph2l+1)
106 : :
107 : : !arrays
108 : : real(dp),intent(in) :: amu(Crystal%ntypat),d2cart(2,3,mpert,3,mpert)
109 : : real(dp),intent(in) :: phfrq(3*natom),epsinf(3,3)
110 : : real(dp),intent(inout) :: displ(2,3*natom,3*natom)
111 : : real(dp),intent(out) :: dielt_rlx(3,3),fact_oscstr(2,3,3*natom)
112 : :
113 : : !Local variables -------------------------
114 : : !scalars
115 : : integer,parameter :: master=0
116 : : integer :: dieflag,idir1,idir2,ifreq,ii,imode,iphl2,nfreq
117 : : integer :: nprocs,my_rank,ncerr
118 : : real(dp) :: afreq,difffr,eps,q2,ucvol
119 : : character(len=500) :: message
120 : : !arrays
121 : : real(dp) :: qphon(3),refl(3)
122 : 26 : real(dp),allocatable :: frdiel(:,:,:),modez(:,:,:),oscstr(:,:,:,:),dielt_modedecompo(:,:,:)
123 : :
124 : : ! *********************************************************************
125 : :
126 : 26 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
127 : :
128 : 26 : dieflag=anaddb_dtset%dieflag
129 : 26 : nfreq=anaddb_dtset%nfreq
130 : :
131 : 26 : ucvol = Crystal%ucvol
132 : :
133 : :
134 : : ! frdiel(3,3,nfreq)= frequency-dependent dielectric tensor
135 : : ! modez=mode effective charge
136 [ - + + - : 78 : ABI_MALLOC(frdiel,(3,3,nfreq))
+ - + - -
+ - + ]
137 [ + - - + : 104 : ABI_MALLOC(oscstr,(2,3,3,3*natom))
+ - + - +
- - + -
+ ]
138 [ - + + - : 78 : ABI_MALLOC(modez,(2,3,3*natom))
+ - + - -
+ - + ]
139 [ - + + - : 78 : ABI_MALLOC(dielt_modedecompo,(3,3,3*natom))
+ - + - -
+ - + ]
140 : :
141 : : !In case only the electronic dielectric constant should be printed
142 [ + + ]: 26 : if (dieflag==2) then
143 : 7 : call ddb_diel_elec(iout,epsinf)
144 : : else
145 : : ! In case the ionic contribution to the dielectric tensor is asked
146 [ + + ]: 19 : if (dieflag/=2 .and. nph2l==0) then
147 : : ! Check if the alignement of phonon modes eigenvector is requested from the input flag alphon;
148 : : ! useful in case of degenerate modes
149 [ + + ]: 10 : if (anaddb_dtset%alphon > 0) then
150 : : write(message, '(3a)' )&
151 : 3 : ' The alphon input variable is non-zero, will mix the degenerate phonon modes',ch10,&
152 : 6 : ' in order to align the mode effective charges with the cartesian axes.'
153 : 3 : call wrtout(std_out,message,'COLL')
154 : 3 : call wrtout(iout,message,'COLL')
155 : 3 : call alignph(amu,displ,d2cart,mpert,natom,Crystal%ntypat,phfrq,Crystal%typat)
156 : : end if ! alignment of the phonon eigenvectors
157 : :
158 : : ! Compute the mode effective charge and oscillator strength
159 : 10 : call ddb_oscstr(displ,d2cart,fact_oscstr,oscstr,modez,iout,mpert,natom,phfrq,ncid,my_rank)
160 : :
161 : : ! Calculation of epsilon_r (Eq.55 PRB 55, 10355)
162 : : ! Check the acousticity of the three lowest modes, assuming
163 : : ! that they are ordered correctly
164 [ - + ]: 10 : if (abs(phfrq(1))>abs(phfrq(4)))then
165 : : ! This means that there is at least one mode with truly negative frequency
166 : : write(message, '(12a,4es16.8)' )&
167 : 0 : 'The lowest mode appears to be a "true" negative mode,',ch10,&
168 : 0 : 'and not an acoustic mode. This precludes the computation',ch10,&
169 : 0 : 'of the frequency-dependent dielectric tensor.',ch10,&
170 : 0 : 'Action : likely there is no action to be taken, although you,',ch10,&
171 : 0 : 'could try to raise your convergence parameters (ecut and k-points).',ch10,&
172 : 0 : 'For your information, here are the four lowest frequencies :',ch10,&
173 : 0 : (phfrq(ii),ii=1,4)
174 : 0 : ABI_ERROR(message)
175 : : end if
176 : :
177 : : ! Compute the relaxed ion dielectric tensor
178 [ + + ]: 40 : do idir1=1,3
179 [ + + ]: 130 : do idir2=1,3
180 : : ! The electronic contribution to epsilon is added
181 : 90 : dielt_rlx(idir1,idir2)=epsinf(idir1,idir2)
182 : : ! calculation of the phonon contribution (ionic) to epsilon
183 [ + + ]: 1227 : do imode=4,3*natom
184 : : ! Note that the acoustic modes are not included : their
185 : : ! oscillator strength should be exactly zero
186 : : ! Also, only the real part of oscstr is taken into account:
187 : : ! the possible imaginary parts of degenerate modes
188 : : ! will cancel.
189 : : dielt_rlx(idir1,idir2)=dielt_rlx(idir1,idir2)+&
190 : 1107 : & oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol
191 : : ! Mode decomposition of epsilon
192 [ - + ]: 1197 : if (dieflag==3)then
193 : 0 : dielt_modedecompo(idir1,idir2,imode)=oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2)*four_pi/ucvol
194 : : endif ! mode decompo of epsilon
195 : : !DEBUG
196 : : ! if(idir1==1 .and. idir2==2)then
197 : : ! write(std_out,'(a,i4,a,3es16.6)')'imode=',imode,' dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)=',&
198 : : !& dielt_rlx(idir1,idir2),oscstr(1,idir1,idir2,imode),phfrq(imode)
199 : : ! endif
200 : : !ENDDEBUG
201 : : end do ! imode
202 : : end do ! idir2
203 : : end do ! idir1
204 : :
205 : : ! Print the electronic dielectric tensor
206 : 10 : call ddb_diel_elec(iout,epsinf)
207 : :
208 : : ! Print the relaxed ion dielectric tensor
209 : 10 : write(message,'(a,a)') ch10,' Relaxed ion dielectric tensor'
210 : 10 : call wrtout(std_out,message,'COLL')
211 : 10 : call wrtout(iout,message,'COLL')
212 : :
213 [ + + ]: 40 : do idir1=1,3
214 : 30 : write(message,'(3f16.8)')(dielt_rlx(idir1,idir2),idir2=1,3)
215 : 30 : call wrtout(std_out,message,'COLL')
216 : 40 : call wrtout(iout,message,'COLL')
217 : : end do
218 : 10 : call wrtout(iout, " ",'COLL')
219 : 10 : call wrtout(std_out, " ",'COLL')
220 : :
221 : : ! Mode decompo of epsilon
222 [ - + ]: 10 : if (dieflag==3) then
223 : 0 : write(message,'(a,a,a,a)') ch10,' Mode by mode decomposition of the ionic dielectric tensor',&
224 : 0 : ch10,' (the electronic contribution is not included)'
225 : 0 : call wrtout(std_out,message,'COLL')
226 : 0 : call wrtout(iout,message,'COLL')
227 [ # # ]: 0 : do imode=4,3*natom
228 : 0 : write(message,'(a,a,i4,a,es14.6,a,a,3f8.3)') ch10,' Mode number ',imode, ' freq = ',phfrq(imode),' Ha', &
229 [ # # # # ]: 0 : ' Mode Z* (|x|, |y|, |z|)', (abs(modez(1,idir1,imode)),idir1=1,3)
230 : 0 : call wrtout(std_out,message,'COLL')
231 : 0 : call wrtout(iout,message,'COLL')
232 [ # # ]: 0 : do idir1=1,3
233 : : ! write(message,'(a,a,i4)') ch10,' Mode number 2',imode
234 : 0 : write(message,'(3f16.8)')(dielt_modedecompo(idir1,idir2,imode),idir2=1,3)
235 : 0 : call wrtout(std_out,message,'COLL')
236 : 0 : call wrtout(iout,message,'COLL')
237 : : end do
238 : : end do
239 : : endif ! mode decompo of epsilon
240 : :
241 : : ! write the relaxed ion dielectric tensor to the netcdf
242 [ + - ]: 10 : if (ncid /= nctk_noid) then
243 : : ncerr = nctk_def_arrays(ncid, [nctkarr_t("emacro_cart_rlx", "dp", &
244 [ + + ]: 20 : "number_of_cartesian_directions, number_of_cartesian_directions")],defmode=.True.)
245 [ - + ]: 10 : NCF_CHECK(ncerr)
246 [ - + ]: 10 : NCF_CHECK(nctk_set_datamode(ncid))
247 [ - + ]: 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "emacro_cart_rlx"), dielt_rlx))
248 : : end if
249 : :
250 : : ! Frequency-dependent dielectric tensor case
251 [ + + ]: 10 : if(dieflag==1) then
252 : 8 : write(message,'(4a)') ch10,ch10,' Calculate the freq. dep. dielectric constant',ch10
253 : 8 : call wrtout(std_out,message,'COLL')
254 : 8 : write(message,'(3a)') ch10,' Frequency dependent dielectric constant:',ch10
255 : 8 : call wrtout(iout,message,'COLL')
256 : :
257 : : ! Check the possibility of asking the frequency-dependent
258 : : ! dielectric tensor (there should be more than one atom in the unit cell)
259 : : ! EB: this check is not in any automatic test
260 [ - + ]: 8 : if(natom==1)then
261 : : write(message, '(6a)' )&
262 : 0 : ' ddb_diel : WARNING -',ch10,&
263 : 0 : ' When there is only one atom in the unit cell',ch10,&
264 : 0 : ' cell, the dielectric tensor is frequency-independent.'!,&
265 : : ! ' Consequently, dieflag has been reset to 2 . '
266 : 0 : call wrtout(std_out,message,'COLL')
267 : 0 : call wrtout(iout,message,'COLL')
268 : : end if
269 : :
270 : 8 : difffr=zero
271 [ + + ]: 8 : if(nfreq>1)difffr=(anaddb_dtset%frmax-anaddb_dtset%frmin)/(nfreq-1)
272 : :
273 [ - + - + ]: 2 : if (nfreq>10 .and. my_rank == master) then
274 : : write(iout, '(a,a,a,a,a,a,a,a)' )&
275 : 2 : ' ddb_diel : the number of frequencies is larger',&
276 : 2 : ' than 10 => I will consider only',ch10,&
277 : 2 : ' the three principal directions, assume that the tensor',ch10,&
278 : 2 : ' is diagonalized, and give dielectric constant and ',ch10,&
279 : 4 : ' reflectivities.'
280 : : write(iout, '(a,a)' )&
281 : 2 : ' Frequency(Hartree) Dielectric constant ',&
282 : 4 : ' Reflectivity '
283 : : write(iout, '(a,a)' )&
284 : 2 : ' x y z',&
285 : 4 : ' x y z'
286 : : end if
287 : :
288 : : ! Loop on frequencies
289 [ + + ]: 214 : do ifreq=1,nfreq
290 : 206 : afreq=anaddb_dtset%frmin+difffr*(ifreq-1)
291 [ + + ]: 824 : do idir1=1,3
292 [ + + ]: 2678 : do idir2=1,3
293 : 1854 : frdiel(idir1,idir2,ifreq)=epsinf(idir1,idir2)
294 [ + + ]: 8817 : do imode=4,3*natom
295 : : ! Note that the acoustic modes are not included : their
296 : : ! oscillator strength should be exactly zero
297 : : ! Also, only the real part of oscstr is taken into account:
298 : : ! the possible imaginary parts of degenerate modes
299 : : ! will cancel.
300 : : frdiel(idir1,idir2,ifreq)=frdiel(idir1,idir2,ifreq)+&
301 : 8199 : & oscstr(1,idir1,idir2,imode) / (phfrq(imode)**2-afreq**2)*four_pi/ucvol
302 : : end do
303 : : end do
304 : : end do
305 : :
306 : : ! Write all this information (actually, there should be a choice of units for the frequencies ...
307 [ + + ]: 214 : if (nfreq>10) then
308 [ + + ]: 800 : do idir1=1,3
309 [ + + ]: 800 : if(frdiel(idir1,idir1,ifreq)<=zero)then
310 : 24 : refl(idir1)=one
311 : : else
312 : : ! See Gervais and Piriou PRB11,3944(1975) [[cite:Gervais1975]].
313 : 576 : refl(idir1)=( (sqrt(frdiel(idir1,idir1,ifreq)) -one) /(sqrt(frdiel(idir1,idir1,ifreq)) +one) )**2
314 : : end if
315 : : end do
316 [ + - ]: 200 : if (my_rank == master) then
317 [ + - + + ]: 800 : write(iout, '(7es12.4)' )afreq,(frdiel(idir1,idir1,ifreq),idir1=1,3),(refl(idir1),idir1=1,3)
318 : : end if
319 : :
320 : : else
321 [ + - ]: 6 : if (my_rank == master) then
322 : 6 : write(iout, '(a,es12.4,a)' )' Full dielectric tensor at frequency',afreq,' Hartree'
323 [ + + ]: 24 : do idir1=1,3
324 : 24 : write(iout, '(3f16.8)' ) (frdiel(idir1,idir2,ifreq),idir2=1,3)
325 : : end do
326 : 6 : write(iout, '(a)' )' '
327 : : end if
328 : : end if
329 : :
330 : : end do ! End of the loop on frequencies
331 : : end if ! End the condition on frequency-dependent dielectric tensor (dieflag=1)
332 : : end if ! dieflag/=2 .and. nph2l==0
333 : : end if
334 : :
335 : : !Calculation of the Lyddane-Sachs-Teller value of the dielectric constant at zero frequency
336 [ + + + - ]: 26 : if(nph2l/=0 .and. dieflag/=2) then
337 : :
338 : : ! Prepare the output
339 : 9 : write(message, '(a,a,a,a)' ) ch10,&
340 : 9 : ' Generalized Lyddane-Sachs-Teller relation at zero frequency :',ch10,&
341 : 18 : ' Direction Dielectric constant'
342 : 9 : call wrtout(std_out,message,'COLL')
343 : 9 : call wrtout(iout,message,'COLL')
344 : :
345 : : ! Examine every wavevector in the phonon list
346 [ + + ]: 23 : do iphl2=1,anaddb_dtset%nph2l
347 : :
348 [ + + ]: 56 : qphon(1:3)=anaddb_dtset%qph2l(1:3,iphl2)
349 : :
350 [ + + + + : 23 : if(abs(qphon(1))>DDB_QTOL .or. abs(qphon(2))>DDB_QTOL .or. abs(qphon(3))>DDB_QTOL)then
+ - ]
351 : 14 : q2=qphon(1)**2+qphon(2)**2+qphon(3)**2
352 : : eps=qphon(1)**2*epsinf(1,1)+qphon(2)**2*epsinf(2,2)+&
353 : : & qphon(3)**2*epsinf(3,3)+ 2* ( qphon(1)*qphon(2)*epsinf(1,2)+&
354 : 14 : & qphon(1)*qphon(3)*epsinf(1,3)+qphon(2)*qphon(3)*epsinf(2,3))
355 : 14 : eps=eps/q2*exp(lst(iphl2)-lst(anaddb_dtset%nph2l+1))
356 [ + - ]: 14 : if (my_rank == master) then
357 : 14 : write(iout, '(3f10.5,f16.8)' )qphon,eps
358 : 14 : write(std_out,'(3f10.5,f16.8)' )qphon,eps
359 : : end if
360 : : end if
361 : : end do
362 : : end if ! End of the condition of nph2l does not vanish for Lyddane-Sachs-Teller
363 : :
364 : :
365 [ - + ]: 26 : ABI_FREE(frdiel)
366 [ - + ]: 26 : ABI_FREE(modez)
367 [ - + ]: 26 : ABI_FREE(oscstr)
368 [ - + ]: 26 : ABI_FREE(dielt_modedecompo)
369 : :
370 : 26 : end subroutine ddb_diel
371 : : !!***
372 : :
373 : :
374 : :
375 : : !!****f* ABINIT/ddb_diel_elec
376 : : !!
377 : : !! NAME
378 : : !! ddb_diel_elec
379 : : !!
380 : : !! FUNCTION
381 : : !! Print the electronic dielectric constant (clamped ions)
382 : : !!
383 : : !! INPUTS
384 : : !! iout=unit number for outputs
385 : : !! epsinf(3,3)= epsilon^infty = electronic contribution to the
386 : : !! dielectric tensor
387 : : !!
388 : : !! OUTPUT
389 : : !!
390 : : !! SOURCE
391 : :
392 : 17 : subroutine ddb_diel_elec(iout,epsinf)
393 : :
394 : : !Arguments -------------------------------
395 : : !scalars
396 : : integer,intent(in) :: iout
397 : : !arrays
398 : : real(dp),intent(in) :: epsinf(3,3)
399 : :
400 : : !Local variables -------------------------
401 : : !scalars
402 : : integer :: idir1,idir2
403 : : character(len=500) :: message
404 : : !arrays
405 : :
406 : 17 : write(message, '(a,a)' ) ch10,' Electronic dielectric tensor'
407 : 17 : call wrtout(std_out,message,'COLL')
408 : 17 : call wrtout(iout,message,'COLL')
409 : :
410 : : !Compute the electronic contribution to the dielectric tensor
411 : : !Needs only the perturbations with E-field from the DDB
412 [ + + ]: 68 : do idir1=1,3
413 : : ! do idir2=1,3
414 : : ! epsinf(idir1,idir2)=d2cart(1,idir1,natom+2,idir2,natom+2)
415 : : ! end do
416 : 51 : write(message, '(3f16.8)' )(epsinf(idir1,idir2),idir2=1,3)
417 : 51 : call wrtout(std_out,message,'COLL')
418 : 68 : call wrtout(iout,message,'COLL')
419 : : end do
420 : 17 : call wrtout(iout, " ",'COLL')
421 : 17 : call wrtout(std_out, " ",'COLL')
422 : :
423 : 17 : end subroutine ddb_diel_elec
424 : : !!***
425 : :
426 : :
427 : :
428 : : !!****f* ABINIT/ddb_oscstr
429 : : !!
430 : : !! NAME
431 : : !! ddb_oscstr
432 : : !!
433 : : !! FUNCTION
434 : : !! Compute the oscillator strength and the mode effective charge
435 : : !!
436 : : !! INPUTS
437 : : !! displ(2,3*natom,3*natom)=
438 : : !! the displacements of atoms in cartesian coordinates.
439 : : !! The first index means either the real or the imaginary part,
440 : : !! The second index runs on the direction and the atoms displaced
441 : : !! The third index runs on the modes.
442 : : !! d2cart(2,3,mpert,3,mpert)=
443 : : !! dynamical matrix, effective charges, dielectric tensor,....
444 : : !! all in cartesian coordinates
445 : : !! iout=unit number for outputs
446 : : !! mpert =maximum number of ipert
447 : : !! natom=number of atoms in unit cell
448 : : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
449 : : !! matrix eigenvalues, except if these are negative, and in this
450 : : !! case, give minus the square root of the absolute value
451 : : !! of the matrix eigenvalues). Hartree units.
452 : : !!
453 : : !! OUTPUT
454 : : !! fact_oscstr(2,3,3*natom)=oscillator strengths for the different eigenmodes,
455 : : !! for different direction of the electric field.
456 : : !! modez(2,3,3*natom)=mode effective charges for the different eigenmodes,
457 : : !! for different directions of the electric field, following
458 : : !! the definition Eq.(53) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
459 : : !! oscstr(2,3,3,3*natom)=oscillator strengths, following
460 : : !! the definition Eq.(54) in PRB55, 10355 (1997) [[cite:Gonze1997a]]
461 : : !!
462 : : !! SOURCE
463 : :
464 : 10 : subroutine ddb_oscstr(displ,d2cart,fact_oscstr,oscstr,modez,iout,mpert,natom,phfrq,ncid,my_rank)
465 : :
466 : : !Arguments -------------------------------
467 : : !scalars
468 : : integer,intent(in) :: iout,mpert,natom,ncid,my_rank
469 : : !arrays
470 : : real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert)
471 : : real(dp),intent(in) :: phfrq(3*natom)
472 : : real(dp),intent(inout) :: displ(2,3*natom,3*natom)
473 : : real(dp),intent(out) :: fact_oscstr(2,3,3*natom),oscstr(2,3,3,3*natom),modez(2,3,3*natom)
474 : :
475 : : !Local variables -------------------------
476 : : !scalars
477 : : integer,parameter :: master=0
478 : : integer :: i1,idir1,idir2,imode,ipert1
479 : : integer :: ncerr
480 : : real(dp) :: usquare
481 : : logical :: t_degenerate
482 : : !arrays
483 : 10 : character(len=1),allocatable :: metacharacter(:)
484 : :
485 : : ! *********************************************************************
486 : :
487 : : ! Get the factors of the oscillator strength, and the mode effective charge for each mode
488 [ + + ]: 163 : do imode=1,3*natom
489 : : usquare=zero
490 [ + + ]: 3438 : do i1=1,3*natom
491 : 3438 : usquare=usquare+displ(1,i1,imode)*displ(1,i1,imode)+displ(2,i1,imode)*displ(2,i1,imode)
492 : : end do
493 [ + + ]: 622 : do idir2=1,3
494 [ + + ]: 1377 : fact_oscstr(:,idir2,imode)=zero
495 [ + + ]: 1377 : modez(:,idir2,imode)=zero
496 [ + + ]: 1989 : do idir1=1,3
497 [ + + ]: 11691 : do ipert1=1,natom
498 : 9855 : i1=idir1+(ipert1-1)*3
499 : : fact_oscstr(:,idir2,imode)=fact_oscstr(:,idir2,imode)+&
500 [ + + ]: 29565 : & displ(:,i1,imode)*d2cart(1,idir1,ipert1,idir2,natom+2)
501 : : modez(:,idir2,imode)=modez(:,idir2,imode)+&
502 : : & displ(:,i1,imode)*&
503 [ + + ]: 30942 : & d2cart(1,idir1,ipert1,idir2,natom+2)/sqrt(usquare)
504 : : end do
505 : : end do
506 : : end do
507 : : end do
508 : :
509 : : ! Examine the degeneracy of each mode. The portability of the echo of the mode effective
510 : : ! charges and oscillator strengths is very hard to guarantee. On the contrary,
511 : : ! the scalar reductions of these quantities are OK.
512 [ + - - + : 40 : ABI_MALLOC(metacharacter,(3*natom))
+ - + - -
+ - + ]
513 [ + + ]: 163 : do imode=1,3*natom
514 : : ! The degenerate modes are not portable
515 : 153 : t_degenerate=.false.
516 [ + + ]: 153 : if(imode>1)then
517 [ + + ]: 143 : if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true.
518 : : end if
519 [ + + ]: 153 : if(imode<3*natom)then
520 [ + + ]: 143 : if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true.
521 : : end if
522 : 153 : metacharacter(imode)=';'
523 [ + + ]: 163 : if(t_degenerate)metacharacter(imode)='-'
524 : : end do
525 : :
526 [ + - ]: 10 : if (my_rank == master) then
527 : : ! Write the mode effective charge for each mode
528 : 10 : write(iout, '(a)' )' '
529 : 10 : write(iout, '(a)' )' Mode effective charges '
530 : 10 : write(iout, '(a)' )' Mode number. x y z length '
531 [ + + ]: 163 : do imode=1,3*natom
532 : 153 : write(iout, '(a,i6,a,4f13.6)' )metacharacter(imode),imode,' ',(modez(1,idir1,imode),idir1=1,3),&
533 : 316 : & (sqrt(modez(1,1,imode)**2+modez(1,2,imode)**2+modez(1,3,imode)**2))
534 : : end do
535 : : end if ! master
536 : :
537 : : ! Get the oscillator strengths
538 [ + + ]: 163 : do imode=1,3*natom
539 [ + + ]: 622 : do idir1=1,3
540 [ + + ]: 1989 : do idir2=1,3
541 : : oscstr(1,idir1,idir2,imode)= &
542 : : & fact_oscstr(1,idir1,imode)*fact_oscstr(1,idir2,imode) +&
543 : 1377 : & fact_oscstr(2,idir1,imode)*fact_oscstr(2,idir2,imode)
544 [ + + ]: 1377 : if(abs(oscstr(1,idir1,idir2,imode))<tol14)oscstr(1,idir1,idir2,imode)=zero
545 : :
546 : : !DEBUG
547 : : ! if(idir1==1 .and. idir2==2)then
548 : : ! write(std_out,'(a,i4,a,5es16.6)')'imode=',imode,&
549 : : !& ' oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)=',&
550 : : !& oscstr(1,idir1,idir2,imode), fact_oscstr(:,idir1,imode),fact_oscstr(:,idir2,imode)
551 : : ! endif
552 : : !ENDDEBUG
553 : :
554 : : oscstr(2,idir1,idir2,imode)= &
555 : : & fact_oscstr(1,idir1,imode)*fact_oscstr(2,idir2,imode) -&
556 : 1377 : & fact_oscstr(2,idir1,imode)*fact_oscstr(1,idir2,imode)
557 [ + - ]: 1836 : if(abs(oscstr(2,idir1,idir2,imode))<tol14)oscstr(2,idir1,idir2,imode)=zero
558 : : end do
559 : : end do
560 : : end do
561 : :
562 [ + - ]: 10 : if (my_rank == master) then
563 : : ! Write the oscillator strength for each mode
564 : 10 : write(iout, '(a)' )' '
565 : 10 : write(iout, '(a)' )' Oscillator strengths (in a.u. ; 1 a.u.=253.2638413 m3/s2). Set to zero if abs()<tol14.'
566 : 10 : write(iout, '(a)' )' Mode number. xx yy zz xy xz yz trace'
567 [ + + ]: 163 : do imode=1,3*natom
568 : : write(iout, '(a,i4,a,7es12.4)' )&
569 [ + + ]: 765 : & metacharacter(imode),imode,' Real ',(oscstr(1,idir1,idir1,imode),idir1=1,3),&
570 : 153 : & oscstr(1,1,2,imode), oscstr(1,1,3,imode),oscstr(1,2,3,imode),&
571 [ + - ]: 918 : & ((oscstr(1,1,1,imode)+oscstr(1,2,2,imode)+oscstr(1,3,3,imode)))
572 : : write(iout, '(a,a,6es12.4)' )&
573 [ + + ]: 765 : & metacharacter(imode),' Imag ',(oscstr(2,idir1,idir1,imode),idir1=1,3),&
574 [ + - ]: 928 : & oscstr(2,1,2,imode),oscstr(2,1,3,imode),oscstr(2,2,3,imode)
575 : : end do
576 : :
577 : : ! write the oscillator strength to the netcdf
578 [ + - ]: 10 : if (ncid /= nctk_noid) then
579 : : ncerr = nctk_def_arrays(ncid, [nctkarr_t("oscillator_strength", "dp", &
580 : : "complex, number_of_cartesian_directions, number_of_cartesian_directions, number_of_phonon_modes")], &
581 [ + + ]: 20 : defmode=.True.)
582 [ - + ]: 10 : NCF_CHECK(ncerr)
583 [ - + ]: 10 : NCF_CHECK(nctk_set_datamode(ncid))
584 [ - + ]: 10 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "oscillator_strength"), oscstr))
585 : : end if
586 : : end if
587 : :
588 [ - + ]: 10 : ABI_FREE(metacharacter)
589 : :
590 : 10 : end subroutine ddb_oscstr
591 : : !!***
592 : :
593 : :
594 : :
595 : : !!****f* ABINIT/alignph
596 : : !!
597 : : !! NAME
598 : : !! alignph
599 : : !!
600 : : !! FUNCTION
601 : : !! Construct linear combinations of the phonon eigendisplacements
602 : : !! of degenerate modes in order to align the mode effective charges
603 : : !! along the axes of the cartesian frame.
604 : : !!
605 : : !! INPUTS
606 : : !! amu(ntypat)=mass of the atoms (atomic mass unit)
607 : : !! displ(2,3*natom,3*natom)=
608 : : !! the displacements of atoms in cartesian coordinates.
609 : : !! The first index means either the real or the imaginary part,
610 : : !! The second index runs on the direction and the atoms displaced
611 : : !! The third index runs on the modes.
612 : : !! d2cart(2,3,mpert,3,mpert)=
613 : : !! dynamical matrix, effective charges, dielectric tensor,....
614 : : !! all in cartesian coordinates
615 : : !! mpert =maximum number of ipert
616 : : !! natom=number of atoms in unit cell
617 : : !! ntypat=number of types of atoms
618 : : !! phfrq(3*natom)=phonon frequencies (square root of the dynamical
619 : : !! matrix eigenvalues, except if these are negative, and in this
620 : : !! case, give minus the square root of the absolute value
621 : : !! of the matrix eigenvalues). Hartree units.
622 : : !! typat(natom)=integer label of each type of atom (1,2,...)
623 : : !!
624 : : !! OUTPUT
625 : : !! displ(2,3*natom,3*natom)=
626 : : !! the displacements of atoms in cartesian coordinates.
627 : : !! The eigendisplacements of degenerate modes have been aligned along
628 : : !! the cartesian axes.
629 : : !!
630 : : !! SOURCE
631 : :
632 : 33 : subroutine alignph(amu,displ,d2cart,mpert,natom,ntypat,phfrq,typat, &
633 : : silent) !optional
634 : :
635 : : !Arguments -------------------------------
636 : : !scalars
637 : : integer,intent(in) :: mpert,natom,ntypat
638 : : integer,optional,intent(in) :: silent
639 : : !arrays
640 : : integer,intent(in) :: typat(natom)
641 : : real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),phfrq(3*natom)
642 : : real(dp),intent(inout) :: displ(2,3*natom,3*natom)
643 : :
644 : : !Local variables -------------------------
645 : : !scalars
646 : : integer,parameter :: master=0
647 : : integer :: i1,idir1,idir2,ii,imode,imodex,imodey,imodez,ipert1,silent_
648 : : real(dp) :: theta
649 : : !arrays
650 : 33 : integer,allocatable :: deg(:)
651 : : real(dp) :: zvec(3,3),zvect(3,3)
652 : 33 : real(dp),allocatable :: modez(:,:,:),modezabs(:),oscstr(:,:,:),vec(:,:),vect(:,:)
653 : :
654 : : ! *********************************************************************
655 : :
656 : : !Reduce verbosity
657 : 33 : silent_=0
658 [ + + ]: 33 : if (present(silent)) silent_=silent
659 : :
660 : : !Get the oscillator strength and mode effective charge for each mode
661 [ + - - + : 132 : ABI_MALLOC(oscstr,(2,3,3*natom))
+ - + - +
- - + -
+ ]
662 [ - + + - : 66 : ABI_MALLOC(modez,(2,3,3*natom))
- + - + ]
663 [ - + + - : 99 : ABI_MALLOC(modezabs,(3*natom))
+ - + - -
+ - + ]
664 [ - + + - : 132 : ABI_MALLOC(vec,(3*natom,3))
+ - + - +
- - + -
+ ]
665 [ - + + - : 66 : ABI_MALLOC(vect,(3*natom,3))
- + - + ]
666 [ - + + - : 99 : ABI_MALLOC(deg,(3*natom))
+ - + - -
+ - + ]
667 : :
668 [ + + ]: 33 : if (silent_/=1) then
669 : 3 : write(std_out,'(a,a)')ch10,' alignph : before modifying the eigenvectors, mode number and mode effective charges :'
670 : : end if
671 [ + + ]: 591 : do imode=1,3*natom
672 : 558 : modezabs(imode)=zero
673 [ + + ]: 1674 : do ii=1,2
674 [ + + ]: 5022 : do idir2=1,3
675 : 3348 : oscstr(ii,idir2,imode)=zero
676 : 3348 : modez(ii,idir2,imode)=zero
677 [ + + ]: 13392 : do idir1=1,3
678 [ + + ]: 72360 : do ipert1=1,natom
679 : 58968 : i1=idir1+(ipert1-1)*3
680 : : oscstr(ii,idir2,imode)=oscstr(ii,idir2,imode)+&
681 : : & displ(ii,i1,imode)*&
682 : 58968 : & d2cart(1,idir1,ipert1,idir2,natom+2)
683 : : modez(ii,idir2,imode)=modez(ii,idir2,imode)+&
684 : : & displ(ii,i1,imode)*&
685 : : & d2cart(1,idir1,ipert1,idir2,natom+2)*&
686 : 69012 : & sqrt(amu(typat(ipert1))*amu_emass)
687 : : end do
688 : : end do
689 [ + + ]: 4464 : if(abs(modez(ii,idir2,imode))>modezabs(imode))modezabs(imode)=abs(modez(ii,idir2,imode))
690 : : end do
691 : : end do
692 [ + + ]: 591 : if (silent_/=1) write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode)
693 : : end do
694 : :
695 : : !Find degenerate modes with non-zero mode effective charge
696 : 33 : imode = 0
697 [ + + ]: 482 : do while (imode < 3*natom)
698 : 449 : imode = imode + 1
699 [ + + ]: 482 : if (imode == 3*natom) then
700 : 10 : deg(imode) = 1
701 [ + + + + : 439 : else if (abs(phfrq(imode) - phfrq(imode+1)) > tol6 .or. modezabs(imode)<tol8 .or. modezabs(imode+1)<tol8) then
- + ]
702 : : ! Differ by phonon frequency or zero mode effective charge
703 : 336 : deg(imode) = 1
704 : : else
705 : 103 : deg(imode) = 2
706 [ + + ]: 103 : if (imode < 3*natom - 1) then
707 [ + + + - ]: 83 : if (abs(phfrq(imode) - phfrq(imode+2)) < tol6 .and. modezabs(imode+2)>tol8) then
708 : 6 : deg(imode) = 3
709 : 6 : imode = imode + 1
710 : : end if
711 : : end if
712 : 103 : imode = imode + 1
713 : : end if
714 : : end do
715 : :
716 : :
717 : : !In case of a degenerate mode, with non-zero mode effective charge, align the mode effective charge vector along
718 : : !the axes of the cartesian frame
719 : 33 : imode = 1
720 [ + + ]: 482 : do while (imode <= 3*natom)
721 : :
722 [ + + ]: 449 : if (silent_/=1) then
723 : 6 : write(std_out,'(a,a,i4,a,i2)')ch10,' Mode number ',imode,' has degeneracy ',deg(imode)
724 : 6 : write(std_out,'(a,3es16.6)') ' Mode effective charge of this mode =',modez(1,:,imode)
725 : : end if
726 : :
727 [ + + ]: 449 : if (deg(imode) == 2) then
728 : :
729 : : ! Optimize on the x direction
730 [ - + ]: 97 : if (silent_/=1) then
731 : 0 : write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1)
732 : : end if
733 [ + + ]: 97 : if (abs(modez(1,1,imode)) > tol8) then
734 : 95 : theta = atan(-modez(1,1,imode+1)/modez(1,1,imode))
735 [ + + ]: 1805 : vec(:,1) = displ(1,:,imode)
736 [ + + ]: 1805 : vec(:,2) = displ(1,:,imode+1)
737 [ + + ]: 1805 : displ(1,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2)
738 [ + + ]: 1805 : displ(1,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2)
739 : :
740 : : ! MR: Rotate also the imaginary part
741 [ + + ]: 1805 : vec(:,1) = displ(2,:,imode)
742 [ + + ]: 1805 : vec(:,2) = displ(2,:,imode+1)
743 [ + + ]: 1805 : displ(2,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2)
744 [ + + ]: 1805 : displ(2,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2)
745 : :
746 : : end if
747 : :
748 [ + + ]: 352 : else if (deg(imode) == 3) then
749 : :
750 [ + - ]: 6 : if (silent_/=1) then
751 : 6 : write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1)
752 : 6 : write(std_out,'(a,3es16.6)') ' Mode effective charge of next-next mode =',modez(1,:,imode+2)
753 : : end if
754 : :
755 : : ! Before mixing them, select the mode-effective charge vectors as being predominently "x", "y" or "z" type.
756 [ # # ]: 0 : if(abs(modez(1,1,imode))>abs(modez(1,2,imode))-tol12 .and. &
757 [ - + ]: 6 : & abs(modez(1,1,imode))>abs(modez(1,3,imode))-tol12) then
758 : 0 : imodex=imode
759 [ # # ]: 0 : if(abs(modez(1,2,imode+1))>abs(modez(1,3,imode+1))-tol12)then
760 : 0 : imodey=imode+1 ; imodez=imode+2
761 : : else
762 : 0 : imodez=imode+1 ; imodey=imode+2
763 : : end if
764 [ - + ]: 6 : else if(abs(modez(1,2,imode))>abs(modez(1,1,imode))-tol12 .and. &
765 [ + - ]: 6 : & abs(modez(1,2,imode))>abs(modez(1,3,imode))-tol12) then
766 : 0 : imodey=imode
767 [ # # ]: 0 : if(abs(modez(1,1,imode+1))>abs(modez(1,3,imode+1))-tol12)then
768 : 0 : imodex=imode+1 ; imodez=imode+2
769 : : else
770 : 0 : imodez=imode+1 ; imodex=imode+2
771 : : end if
772 : : else
773 : 6 : imodez=imode
774 [ - + ]: 6 : if(abs(modez(1,1,imode+1))>abs(modez(1,2,imode+1))-tol12)then
775 : 0 : imodex=imode+1 ; imodey=imode+2
776 : : else
777 : 6 : imodey=imode+1 ; imodex=imode+2
778 : : end if
779 : : end if
780 [ + + ]: 42 : vec(:,1)=displ(1,:,imodex)
781 [ + + ]: 42 : vec(:,2)=displ(1,:,imodey)
782 [ + + ]: 42 : vec(:,3)=displ(1,:,imodez)
783 [ + + ]: 24 : zvec(:,1)=modez(1,:,imodex)
784 [ + + ]: 24 : zvec(:,2)=modez(1,:,imodey)
785 [ + + ]: 24 : zvec(:,3)=modez(1,:,imodez)
786 : :
787 : :
788 : : ! Optimize along x : does the first vector has a component along x ?
789 [ + - ]: 6 : if (abs(zvec(1,1)) > tol8) then
790 : : ! Optimize on the (1,2) pair of modes along x
791 : 6 : theta = atan(-zvec(1,2)/zvec(1,1))
792 : 6 : zvect(:,:)=zvec(:,:)
793 [ + + ]: 24 : zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,2)
794 [ + + ]: 24 : zvec(:,2) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,2)
795 [ + + + + ]: 132 : vect(:,:)=vec(:,:)
796 [ + + ]: 42 : vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,2)
797 [ + + ]: 42 : vec(:,2) = sin(theta)*vect(:,1) + cos(theta)*vect(:,2)
798 : : ! Optimize on the (1,3) pair of modes along x
799 : 6 : theta = atan(-zvec(1,3)/zvec(1,1))
800 : 6 : zvect(:,:)=zvec(:,:)
801 [ + + ]: 24 : zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3)
802 [ + + ]: 24 : zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3)
803 [ + + + + ]: 132 : vect(:,:)=vec(:,:)
804 [ + + ]: 42 : vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3)
805 [ + + ]: 42 : vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3)
806 [ + - ]: 6 : if (abs(zvec(2,2)) > tol8) then
807 : : ! Optimize on the (2,3) pair of modes along y
808 : 6 : theta = atan(-zvec(2,3)/zvec(2,2))
809 : 6 : zvect(:,:)=zvec(:,:)
810 [ + + ]: 24 : zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
811 [ + + ]: 24 : zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
812 [ + + + + ]: 132 : vect(:,:)=vec(:,:)
813 [ + + ]: 42 : vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
814 [ + + ]: 42 : vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
815 : : end if
816 : : ! Likely, the remaining is not needed ... because the vectors have been ordered in x, y, and z major component ...
817 : : ! Optimize along x : does the second vector has a component along x ?
818 [ # # ]: 0 : else if(abs(zvec(1,2)) > tol8) then
819 : : ! Optimize on the (2,3) pair of modes along x
820 : 0 : theta = atan(-zvec(1,3)/zvec(1,2))
821 : 0 : zvect(:,:)=zvec(:,:)
822 [ # # ]: 0 : zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
823 [ # # ]: 0 : zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
824 [ # # # # ]: 0 : vect(:,:)=vec(:,:)
825 [ # # ]: 0 : vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
826 [ # # ]: 0 : vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
827 : : ! Optimize on the (1,3) pair of modes along y
828 [ # # ]: 0 : if (abs(zvec(2,1)) > tol8) then
829 : 0 : theta = atan(-zvec(2,3)/zvec(2,1))
830 : 0 : zvect(:,:)=zvec(:,:)
831 [ # # ]: 0 : zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3)
832 [ # # ]: 0 : zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3)
833 [ # # # # ]: 0 : vect(:,:)=vec(:,:)
834 [ # # ]: 0 : vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3)
835 [ # # ]: 0 : vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3)
836 : : end if
837 : : ! We are left with the pair of vectors (2,3)
838 [ # # ]: 0 : else if (abs(zvec(2,2)) > tol8) then
839 : : ! Optimize on the (2,3) pair of modes along y
840 : 0 : theta = atan(-zvec(2,3)/zvec(2,2))
841 : 0 : zvect(:,:)=zvec(:,:)
842 [ # # ]: 0 : zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3)
843 [ # # ]: 0 : zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3)
844 [ # # # # ]: 0 : vect(:,:)=vec(:,:)
845 [ # # ]: 0 : vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3)
846 [ # # ]: 0 : vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3)
847 : : end if
848 : :
849 [ + + ]: 42 : displ(1,:,imodex)=vec(:,1)
850 [ + + ]: 42 : displ(1,:,imodey)=vec(:,2)
851 [ + + ]: 42 : displ(1,:,imodez)=vec(:,3)
852 : :
853 : : ! Previous coding, from Marek. Apparently, break the orthogonalization of vectors ...
854 : : ! do ii = 1,3
855 : : ! coeff(:) = 0._dp
856 : : ! if (ii == 1) then
857 : : ! jj = 2 ; kk = 3
858 : : ! else if (ii == 2) then
859 : : ! jj = 1 ; kk = 3
860 : : ! else
861 : : ! jj = 1 ; kk = 2
862 : : ! end if
863 : : ! coeff(ii) = 1._dp
864 : : ! c1 = modez(1,jj,imode+ii-1)
865 : : ! c2 = modez(1,jj,imode+jj-1)
866 : : ! c3 = modez(1,jj,imode+kk-1)
867 : : ! c4 = modez(1,kk,imode+ii-1)
868 : : ! c5 = modez(1,kk,imode+jj-1)
869 : : ! c6 = modez(1,kk,imode+kk-1)
870 : : ! dtm = c2*c6 - c3*c5
871 : : ! if (abs(dtm) > tol8) then
872 : : ! coeff(jj) = (c3*c4 - c1*c6)/dtm
873 : : ! coeff(kk) = (c1*c5 - c2*c4)/dtm
874 : : ! end if
875 : : ! mod_ = sqrt(1._dp + coeff(jj)*coeff(jj) + coeff(kk)*coeff(kk))
876 : : ! coeff(:) = coeff(:)/mod_
877 : : ! displ(1,:,imode+ii-1) = coeff(1)*vec(1,:) + coeff(2)*vec(2,:) + &
878 : : !& coeff(3)*vec(3,:)
879 : : ! end do
880 : :
881 : : end if ! if deg mode
882 : :
883 : 449 : imode = imode + deg(imode)
884 : :
885 : : end do
886 : :
887 [ + + ]: 33 : if (silent_/=1) then
888 : 3 : write(std_out,'(a,a)')ch10,' alignph : after modifying the eigenvectors, mode number and mode effective charges :'
889 : : end if
890 [ + + ]: 591 : do imode=1,3*natom
891 [ + + ]: 1674 : do ii=1,2
892 [ + + ]: 5022 : do idir2=1,3
893 : 3348 : modez(ii,idir2,imode)=zero
894 [ + + ]: 14508 : do idir1=1,3
895 [ + + ]: 72360 : do ipert1=1,natom
896 : 58968 : i1=idir1+(ipert1-1)*3
897 : : modez(ii,idir2,imode)=modez(ii,idir2,imode)+&
898 : : & displ(ii,i1,imode)*&
899 : : & d2cart(1,idir1,ipert1,idir2,natom+2)*&
900 : 69012 : & sqrt(amu(typat(ipert1))*amu_emass)
901 : : end do
902 : : end do
903 : : end do
904 : : end do
905 [ + + ]: 591 : if (silent_/=1) write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode)
906 : : end do
907 : :
908 [ - + ]: 33 : ABI_FREE(deg)
909 [ - + ]: 33 : ABI_FREE(oscstr)
910 [ - + ]: 33 : ABI_FREE(modez)
911 [ - + ]: 33 : ABI_FREE(modezabs)
912 [ - + ]: 33 : ABI_FREE(vec)
913 [ - + ]: 33 : ABI_FREE(vect)
914 : :
915 : 33 : end subroutine alignph
916 : : !!***
917 : :
918 : : end module m_ddb_diel
919 : : !!***
|
