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