Branch data Line data Source code
1 : : !!****m* ABINIT/m_sorth_ph
2 : : !! NAME
3 : : !!
4 : : !! FUNCTION
5 : : !!
6 : : !! COPYRIGHT
7 : : !! Copyright (C) 2008-2026 ABINIT group (MVer, FDortu, MVeithen)
8 : : !! This file is distributed under the terms of the
9 : : !! GNU General Public License, see ~abinit/COPYING
10 : : !! or http://www.gnu.org/copyleft/gpl.txt .
11 : : !!
12 : : !! SOURCE
13 : :
14 : : #if defined HAVE_CONFIG_H
15 : : #include "config.h"
16 : : #endif
17 : :
18 : : #include "abi_common.h"
19 : :
20 : : module m_sortph
21 : :
22 : : use defs_basis
23 : : use m_abicore
24 : : use m_errors
25 : :
26 : : use m_io_tools, only : open_file
27 : :
28 : : implicit none
29 : :
30 : : private
31 : :
32 : : complex(dp),save,allocatable :: eigvecLast(:,:)
33 : :
34 : : public :: end_sortph
35 : : public :: sortph
36 : :
37 : : ! Logical units used to write data.
38 : : integer,private,save :: udispl=-1,ufreq=-1
39 : : !!***
40 : :
41 : : contains
42 : :
43 : : !!****f* m_sortph/end_sortph
44 : : !! NAME
45 : : !! end_sortph
46 : : !!
47 : : !! FUNCTION
48 : : !! Deallocate array for sortph
49 : : !!
50 : : !! INPUTS
51 : : !!
52 : : !! OUTPUT
53 : : !! Only deallocation
54 : : !!
55 : : !! NOTES
56 : : !!
57 : : !! SOURCE
58 : 79 : subroutine end_sortph()
59 : :
60 [ + + ]: 79 : if (allocated(eigvecLast)) then
61 : 12 : ABI_FREE(eigvecLast)
62 : : end if
63 : :
64 [ + + ]: 79 : if (ufreq /= -1) then
65 : 12 : close(ufreq); ufreq = -1
66 : : end if
67 [ + + ]: 79 : if (udispl /= -1) then
68 : 12 : close(udispl); udispl = -1
69 : : end if
70 : :
71 : 79 : end subroutine end_sortph
72 : : !!***
73 : :
74 : : !!****f* m_sortph/sortph
75 : : !! NAME
76 : : !! sortph
77 : : !!
78 : : !! FUNCTION
79 : : !! Sort the energies in order to have fine phonon dispersion curves
80 : : !! It is best not to include the gamma point in the list
81 : : !!
82 : : !! MODIFIED
83 : : !! Takeshi Nishimatsu
84 : : !!
85 : : !! INPUTS
86 : : !! eigvec(2*3*natom*3*natom)= contain
87 : : !! the eigenvectors of the dynamical matrix.
88 : : !! displ(2*3*natom*3*natom)= contain
89 : : !! the displacements of atoms in cartesian coordinates.
90 : : !! The first index means either the real or the imaginary part,
91 : : !! The second index runs on the direction and the atoms displaced
92 : : !! The third index runs on the modes.
93 : : !! filnam=name of output files
94 : : !! hacmm1,hartev,harthz,xkb= different conversion factors
95 : : !! natom= number of atom
96 : : !! phfrq(3*natom)= phonon frequencies in Hartree
97 : : !!
98 : : !! OUTPUT
99 : : !! (only writing ?)
100 : : !!
101 : : !! NOTES
102 : : !! Called by one processor only
103 : : !!
104 : : !! SOURCE
105 : :
106 : 588 : subroutine sortph(eigvec,displ,filnam, natom,phfrq)
107 : :
108 : : !Arguments -----------------------------------
109 : : !scalars
110 : : integer,intent(in) :: natom
111 : : character(len=*),intent(in) :: filnam
112 : : !arrays
113 : : real(dp),intent(in) :: eigvec(2,3,natom,3,natom)
114 : : real(dp),intent(in) :: displ(2*3*natom*3*natom)
115 : : real(dp),intent(in) :: phfrq(3*natom)
116 : :
117 : : !Local variables-------------------------------
118 : : !scalars
119 : : integer :: iatom,imode,j,idir1,idir2,ipert1,ipert2,i1,i2
120 : : character(len=fnlen) :: file_displ,file_freq
121 : : character(len=20) :: fmt_phfrq
122 : : character(len=500) :: msg
123 : : !arrays
124 : 1176 : logical :: mask(3*natom)
125 : 1176 : real(dp) :: phfrqNew(3*natom)
126 : 1176 : complex(dp) :: displIn(3*natom,3*natom)
127 : 1176 : complex(dp) :: displNew(3*natom,3*natom)
128 : 1176 : complex(dp) :: eigvecIn(3*natom,3*natom)
129 : 1176 : complex(dp) :: eigvecNew(3*natom,3*natom)
130 : 1176 : complex(dp) :: transpose_eigvec(3*natom,3*natom)
131 : 1176 : real(dp) :: abs_similarity(3*natom,3*natom) !|<displNew|displLast>|
132 : : ! *********************************************************************
133 : :
134 [ + + ]: 2276 : do ipert2=1,natom
135 [ + + ]: 7340 : do idir2=1,3
136 : 5064 : i2=idir2+(ipert2-1)*3
137 [ + + ]: 23390 : do ipert1=1,natom
138 [ + + ]: 71616 : do idir1=1,3
139 : 49914 : i1=idir1+(ipert1-1)*3
140 : 49914 : eigvecIn(i1,i2)=cmplx(eigvec(1,idir1,ipert1,idir2,ipert2),eigvec(2,idir1,ipert1,idir2,ipert2))
141 : 66552 : displIn(i1,i2)=cmplx(displ(1+2*(i1-1)+2*3*natom*(i2-1)),displ(2+2*(i1-1)+2*3*natom*(i2-1)))
142 : : end do
143 : : end do
144 : : end do
145 : : end do
146 : :
147 [ + + ]: 588 : if(.not.allocated(eigvecLast)) then
148 [ + - ]: 12 : file_freq = trim(filnam)//".freq" !---------------------------------------------------
149 : 12 : write(std_out,'(a,a)' )' sortph : opening file ',trim(file_freq)
150 [ - + ]: 12 : if (open_file(file_freq,msg,newunit=ufreq,STATUS='replace',ACTION='write') /= 0) then
151 : 0 : ABI_ERROR(msg)
152 : : end if
153 [ + - ]: 12 : file_displ = trim(filnam)//".displ" !--------------------------------------------------
154 : 12 : write(std_out,'(a,a)' )' sortph : opening file ',trim(file_displ)
155 [ - + ]: 12 : if (open_file(file_displ,msg,newunit=udispl,STATUS='replace',ACTION='write') /= 0) then
156 : 0 : ABI_ERROR(msg)
157 : : end if
158 [ + - - + : 60 : ABI_MALLOC(eigvecLast,(3*natom,3*natom))
+ - - + +
- + - + -
- + - + -
+ ]
159 [ + + ]: 108 : phfrqNew(:) = phfrq(:)
160 [ + + + + ]: 990 : displNew(:,:) = displIn(:,:)
161 [ + + + + ]: 990 : eigvecNew(:,:) = eigvecIn(:,:)
162 : : else
163 : : ! Avoid gfortran 4.2.1 bug, with which you CANNOT conjg(transpose(displ))
164 [ + + + + ]: 54576 : transpose_eigvec = transpose(eigvecIn)
165 [ - + + - : 1247040 : abs_similarity = abs(matmul(conjg(transpose_eigvec),eigvecLast))
- + + - +
- + - - +
- + + + +
+ + + + +
+ + + - +
- - - - -
- - + + +
+ + - ]
166 [ + + ]: 5544 : mask(:) = .true.
167 [ + + ]: 5544 : phfrqNew(:) = phfrq(:)
168 [ + + + + ]: 54576 : displNew(:,:) = displIn(:,:)
169 [ + + + + ]: 54576 : eigvecNew(:,:) = eigvecIn(:,:)
170 : : end if
171 : :
172 : :
173 : : !Write frequencies in a file
174 : 588 : write(fmt_phfrq,'(a,i3,a)') '(', 3*natom, 'e18.10)'
175 : 588 : write(ufreq,fmt_phfrq) (phfrqNew(j),j=1,3*natom)
176 : :
177 : : !write displacements in a file
178 : : ! NB: sqrt still returns a complex number could be using modulus or something simpler
179 [ + + ]: 5652 : do imode=1,3*natom
180 [ + + ]: 22290 : do iatom=1,natom
181 : : write(udispl,'(e18.10)') &
182 : : real(sqrt(displNew(3*(iatom-1)+1,imode) * conjg(displNew(3*(iatom-1)+1,imode)) + &
183 : : & displNew(3*(iatom-1)+2,imode) * conjg(displNew(3*(iatom-1)+2,imode)) + &
184 : 21702 : & displNew(3*(iatom-1)+3,imode) * conjg(displNew(3*(iatom-1)+3,imode)) ))
185 : : end do
186 : : end do
187 : :
188 [ + + + + ]: 55566 : eigvecLast(:,:) = eigvecNew(:,:)
189 : :
190 : 588 : end subroutine sortph
191 : : !!***
192 : :
193 [ + - ]: 576 : end module m_sortph
194 : : !!***
|
