TABLE OF CONTENTS


ABINIT/tdep [ Programs ]

[ Top ] [ Programs ]

NAME

 tdep

FUNCTION

 Calculations of phonons using molecular dynamic simulations

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (FB,JB)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt

NOTES

  The input files are input.in, xred.dat, fcart.dat and etot.dat
  See the examples in the test directory

TODO

  A lot of things

INPUTS

  (main routine)

OUTPUT

  (main routine)

PARENTS

CHILDREN

      tdep_build_phijnn,tdep_calc_dij,tdep_calc_elastic,tdep_calc_model
      tdep_calc_moorepenrose,tdep_calc_phdos,tdep_calc_phijfcoeff
      tdep_calc_thermo,tdep_destroy_shell,tdep_init_crystal,tdep_init_ddb
      tdep_init_eigen2nd,tdep_init_ifc,tdep_init_shell2at,tdep_make_latt
      tdep_make_qptpath,tdep_make_sym,tdep_matchideal2average
      tdep_print_aknowledgments,tdep_readecho,tdep_write_dij,tdep_write_yaml
      xmpi_end,xmpi_init

SOURCE

 43 #if defined HAVE_CONFIG_H
 44 #include "config.h"
 45 #endif
 46 
 47 #include "abi_common.h"
 48 
 49 program tdep
 50 
 51   use defs_basis
 52   use m_abicore
 53   use m_phonons
 54   use m_xmpi,             only : xmpi_init, xmpi_end
 55   use m_ifc,              only : ifc_type
 56   use m_crystal,          only : crystal_t
 57   use m_ddb,              only : ddb_type
 58   use m_dynmat,           only : ftifc_r2q, ftifc_q2r, asrif9
 59   use m_tdep_abitypes,    only : tdep_init_crystal, tdep_init_ifc, tdep_init_ddb, tdep_write_ifc
 60   use m_tdep_phij,        only : tdep_calc_phijfcoeff, tdep_build_phijNN, tdep_calc_dij, tdep_write_dij, &
 61 &                                Eigen_Variables_type, tdep_init_eigen2nd, tdep_destroy_eigen2nd, tdep_write_yaml
 62   use m_tdep_latt,        only : tdep_make_latt, Lattice_Variables_type
 63   use m_tdep_sym,         only : tdep_make_sym, Symetries_Variables_type
 64   use m_tdep_readwrite,   only : tdep_print_Aknowledgments, tdep_ReadEcho, Input_Variables_type
 65   use m_tdep_utils,       only : Coeff_Moore_type, tdep_calc_MoorePenrose, tdep_MatchIdeal2Average, tdep_calc_model
 66   use m_tdep_qpt,         only : tdep_make_qptpath, Qpoints_type
 67   use m_tdep_phdos,       only : tdep_calc_phdos,tdep_calc_elastic,tdep_calc_thermo
 68   use m_tdep_shell,       only : Shell_Variables_type, tdep_init_shell2at, tdep_init_shell3at, tdep_destroy_shell
 69 #ifdef HAVE_NETCDF
 70   use netcdf
 71 #endif
 72   use m_io_tools
 73 
 74 !This section has been created automatically by the script Abilint (TD).
 75 !Do not modify the following lines by hand.
 76 #undef ABI_FUNC
 77 #define ABI_FUNC 'tdep'
 78 !End of the abilint section
 79 
 80   implicit none
 81 
 82   integer :: natom,natom_unitcell,ntotcoeff,nshell_max
 83   integer :: order,stdout,norder,iqpt
 84   double precision :: U0,DeltaFree_AH2
 85   double precision, allocatable :: ucart(:,:,:),proj(:,:,:),proj_tmp(:,:,:),Forces_TDEP(:)
 86 !FB  double precision, allocatable :: fcoeff(:,:),Phij_coeff(:,:),Forces_MD(:),Phij_NN(:,:)
 87   double precision, allocatable :: Phij_coeff(:,:),Forces_MD(:),Phij_NN(:,:)
 88   double precision, allocatable :: distance(:,:,:),Rlatt_cart(:,:,:),Rlatt4Abi(:,:,:)
 89   double precision, allocatable :: omega (:)
 90   double precision :: qpt_cart(3)
 91   double complex  , allocatable :: dij(:,:),eigenV(:,:)
 92   type(phonon_dos_type) :: PHdos
 93   type(Input_Variables_type) :: InVar
 94   type(Lattice_Variables_type) :: Lattice
 95   type(Symetries_Variables_type) :: Sym
 96   type(Qpoints_type) :: Qpt
 97   type(ifc_type) :: Ifc
 98   type(ddb_type) :: DDB
 99   type(crystal_t) :: Crystal
100   type(Shell_Variables_type) :: Shell2at
101   type(Coeff_Moore_type) :: CoeffMoore
102   type(Eigen_Variables_type) :: Eigen2nd
103 
104 !******************************************************************
105 
106 !==========================================================================================
107 !===================== Initialization & Reading  ==========================================
108 !==========================================================================================
109  call xmpi_init()
110 
111 ! Read input values from the input.in input file
112  call tdep_ReadEcho(InVar)
113 
114 ! Initialize basic quantities  
115  natom         =InVar%natom
116  natom_unitcell=InVar%natom_unitcell
117  stdout        =InVar%stdout  
118  nshell_max    =500
119 
120 !==========================================================================================
121 !============== Define the ideal lattice, symmetries and Brillouin zone ===================
122 !==========================================================================================
123 ! Define all the quantities needed to buid the lattice (rprim*, acell*, brav*...)
124  call tdep_make_latt(InVar,Lattice)
125 
126 ! Compute all the symmetries coming from the bravais lattice
127  call tdep_make_sym(Invar,Lattice,Sym)
128 
129 ! Initialize the Brillouin zone and compute the q-points path 
130  call tdep_make_qptpath(InVar,Lattice,Qpt)
131 
132 !==========================================================================================
133 !======== 1/ Determine ideal positions and distances ======================================
134 !======== 2/ Find the matching between the ideal and average ==============================
135 !========   (from the MD simulations) positions. ==========================================
136 !======== 3/ Find the symmetry operation between the reference and image bonds =============
137 !======== 4/ Write output quantities needed to visualize the neighbouring distances =======
138 !==========================================================================================
139  ABI_MALLOC(Rlatt4Abi ,(3,natom_unitcell,natom))   ; Rlatt4Abi (:,:,:)=0.d0
140  ABI_MALLOC(distance,(natom,natom,4))              ; distance(:,:,:)=0.d0
141  ABI_MALLOC(Rlatt_cart,(3,natom_unitcell,natom))   ; Rlatt_cart(:,:,:)=0.d0
142  ABI_MALLOC(ucart,(3,natom,InVar%nstep))    ; ucart(:,:,:)=0.d0
143  ABI_MALLOC(Forces_MD,(3*natom*InVar%nstep)); Forces_MD(:)=0.d0
144 
145  write(InVar%stdout,*) "Matching structure"
146  call flush_unit(InVar%stdout)
147  call tdep_MatchIdeal2Average(distance,Forces_MD,InVar,Lattice,Rlatt_cart,Rlatt4Abi,Sym,ucart)
148  call flush_unit(InVar%stdout)
149 
150 !==========================================================================================
151 !============== Initialize Crystal and DDB ABINIT Datatypes ===============================
152 !==========================================================================================
153  call tdep_init_crystal(Crystal,InVar,Lattice,Sym)
154  call tdep_init_ddb(DDB,InVar,Lattice)
155 
156 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=
157 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=
158 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=# CALCULATION OF THE 2nd ORDER =#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=
159 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=
160 !#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=
161  order=2
162  norder=3**order
163  write(InVar%stdout,*) ' '
164  write(InVar%stdout,*) '#############################################################################'
165  write(InVar%stdout,*) '############################## SECOND ORDER  ################################'
166  write(InVar%stdout,*) '################ Now, find the number of coefficients for ###################'
167  write(InVar%stdout,*) '########################## a reference interaction ##########################'
168  write(InVar%stdout,*) '#############################################################################'
169  call flush_unit(InVar%stdout)
170  
171 !==========================================================================================
172 !============== Initialize the Shell2at datatype ==========================================
173 !==========================================================================================
174  ABI_MALLOC(proj_tmp,(norder,norder,nshell_max)) ; proj_tmp(:,:,:)=0.d0
175  call tdep_init_shell2at(distance,InVar,norder,nshell_max,ntotcoeff,order,proj_tmp,Shell2at,Sym)
176  ABI_MALLOC(proj    ,(norder,norder,Shell2at%nshell)) ; proj(:,:,:)=0.d0
177  proj = reshape (proj_tmp, (/ norder,norder,Shell2at%nshell /)) 
178  ABI_FREE(proj_tmp)
179 
180 !==========================================================================================
181 !============== Initialize the IFC Abinit datatype ========================================
182 !==========================================================================================
183  ABI_MALLOC(Phij_NN,(3*natom,3*natom)) ; Phij_NN(:,:)=0.d0
184  call tdep_init_ifc(Crystal,DDB,Ifc,InVar,Lattice,Phij_NN,Rlatt4Abi,Shell2at,Sym)
185 
186 !==========================================================================================
187 !============= Build fcoeff, needed for the Moore-Penrose method just below ===============
188 !==========================================================================================
189  ABI_MALLOC(CoeffMoore%fcoeff,(3*natom*InVar%nstep,ntotcoeff)); CoeffMoore%fcoeff(:,:)=0.d0 
190  if (InVar%ReadIFC.ne.1) then
191    call tdep_calc_phijfcoeff(InVar,ntotcoeff,proj,Shell2at,Sym,ucart,CoeffMoore%fcoeff)
192  end if
193 
194 !==========================================================================================
195 !============= Compute the pseudo inverse using the Moore-Penrose method ==================
196 !==========================================================================================
197  ABI_MALLOC(Phij_coeff,(ntotcoeff,1)); Phij_coeff(:,:)=0.d0
198  if (InVar%ReadIFC.ne.1) then
199    call tdep_calc_MoorePenrose(Forces_MD,CoeffMoore,InVar,ntotcoeff,Phij_coeff)
200  end if
201 
202 !==========================================================================================
203 !============= Reorganize the IFC coefficients into the whole Phij_NN matrix ==============
204 !==========================================================================================
205  if (InVar%ReadIFC.ne.1) then
206    call tdep_build_phijNN(distance,InVar,ntotcoeff,proj,Phij_coeff,Phij_NN,Shell2at,Sym)
207  end if
208  ABI_FREE(Phij_coeff)
209  ABI_FREE(proj)
210 
211 !==========================================================================================
212 !===================== Compute the phonons density of states ==============================
213 !==========================================================================================
214  call tdep_calc_phdos(Crystal,Ifc,InVar,Lattice,natom,natom_unitcell,Phij_NN,PHdos,Qpt,Rlatt4Abi,Shell2at,Sym)
215  ABI_FREE(Rlatt4Abi)
216  call tdep_destroy_shell(natom,order,Shell2at)
217 
218 !==========================================================================================
219 !============= Compute the dynamical matrix and phonon spectrum ===========================
220 !================ then print Dij, omega, eigenvectors =====================================
221 !==========================================================================================
222  write(stdout,*)' '
223  write(stdout,*) '#############################################################################'
224  write(stdout,*) '######################## Dynamical matrix ###################################'
225  write(stdout,*) '#############################################################################'
226  open(unit=53,file=trim(InVar%output_prefix)//'omega.dat')
227  open(unit=52,file=trim(InVar%output_prefix)//'dij.dat')
228  open(unit=51,file=trim(InVar%output_prefix)//'eigenvectors.dat')
229  ABI_MALLOC(dij   ,(3*InVar%natom_unitcell,3*InVar%natom_unitcell)) 
230  ABI_MALLOC(eigenV,(3*InVar%natom_unitcell,3*InVar%natom_unitcell)) 
231  ABI_MALLOC(omega,(3*InVar%natom_unitcell))
232  call tdep_init_eigen2nd(Eigen2nd,InVar%natom_unitcell,Qpt%nqpt)
233  do iqpt=1,Qpt%nqpt
234    dij(:,:)=zero ; eigenV(:,:)=zero ; omega(:)=zero
235    qpt_cart(:)=Qpt%qpt_cart(:,iqpt)
236    call tdep_calc_dij (dij,eigenV,iqpt,InVar,Lattice,omega,Phij_NN,qpt_cart,Rlatt_cart)
237    call tdep_write_dij(dij,eigenV,iqpt,InVar,Lattice,omega,qpt_cart)
238    Eigen2nd%eigenval(:,iqpt)=  omega(:)
239    Eigen2nd%eigenvec(:,:,iqpt)=eigenV(:,:)
240  end do  
241  ABI_FREE(dij)
242  ABI_FREE(eigenV)
243  ABI_FREE(omega)
244  close(53)
245  close(52)
246  close(51)
247  call tdep_write_yaml(Eigen2nd,Qpt,InVar%output_prefix)
248  write(InVar%stdout,'(a)') ' See the dij.dat, omega.dat and eigenvectors files'
249 !==========================================================================================
250 !===================== Compute the elastic constants ======================================
251 !==========================================================================================
252  call tdep_calc_elastic(Phij_NN,distance,InVar,Lattice)
253 
254 !==========================================================================================
255 !=========== Compute U_0, the "free energy" and the forces (from the model) ===============
256 !==========================================================================================
257  ABI_MALLOC(Forces_TDEP,(3*InVar%natom*InVar%nstep)); Forces_TDEP(:)=0.d0 
258  call tdep_calc_model(DeltaFree_AH2,distance,Forces_MD,Forces_TDEP,InVar,Phij_NN,ucart,U0) 
259 
260 !==========================================================================================
261 !===================== Compute the thermodynamical quantities =============================
262 !==========================================================================================
263  call tdep_calc_thermo(DeltaFree_AH2,InVar,PHdos,U0)
264 
265  !if (InVar%Order==2) then
266  !  call tdep_print_Aknowledgments(InVar)
267  !  call xmpi_end()
268  !  stop
269  !end if  
270 
271  ABI_FREE(Rlatt_cart)
272  ABI_FREE(distance)
273  ABI_FREE(Forces_MD)
274  ABI_FREE(Forces_TDEP)
275  ABI_FREE(Phij_NN)
276  ABI_FREE(ucart)
277 !==========================================================================================
278 !================= Write the last informations (aknowledgments...)  =======================
279 !==========================================================================================
280  call tdep_print_Aknowledgments(InVar)
281  call xmpi_end()
282 
283  end program tdep