Line data Source code
1 : !!****m* ABINIT/m_tdep_model
2 : !! NAME
3 : !! m_tdep_model
4 : !!
5 : !! FUNCTION
6 : !! This module contains the TDEP Model data type
7 : !! which holds the IFC at all orders and derived quantities.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_tdep_model
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_xmpi
29 : use m_io_tools
30 : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
31 : use m_tdep_sampling, only : tdep_Sampling_type
32 : use m_tdep_shell, only : Shell_type
33 :
34 : implicit none
35 :
36 : !----------------------------------------------------------------------
37 :
38 : !!***
39 : !!****t* m_tdep_model/Phi2_type
40 : !! NAME
41 : !! Phi2_type
42 : !!
43 : !! FUNCTION
44 : !! Second-order IFC
45 : !!
46 : !! SOURCE
47 :
48 : type Phi2_type
49 :
50 : double precision, allocatable :: SR(:,:)
51 : ! SR(3*natom, 3*natom)
52 : ! Short-range part of the second-order IFC
53 :
54 : double precision, allocatable :: LR(:,:)
55 : ! LR(3*natom, 3*natom)
56 : ! Long-range part of the second-order IFC
57 :
58 : double precision, allocatable :: Tot(:,:)
59 : ! Tot(3*natom, 3*natom)
60 : ! Second-order IFC (SR + LR)
61 :
62 :
63 : end type Phi2_type
64 :
65 :
66 : !----------------------------------------------------------------------
67 :
68 : !!***
69 : !!****t* m_tdep_model/tdep_Model_type
70 : !! NAME
71 : !! tdep_Model_type
72 : !!
73 : !! FUNCTION
74 : !! Holds the result of the TDEP fit, that is, the IFCs and resulting forces
75 : !! and energy contributions.
76 : !!
77 : !! SOURCE
78 :
79 : type,public :: tdep_Model_type
80 :
81 : integer :: order
82 : ! Maximum IFC order (2, 3, or 4)
83 :
84 : integer :: natom
85 : ! Number of atoms
86 :
87 : integer :: natom_unitcell
88 : ! Number of atoms in the unitcell
89 :
90 : integer :: my_nstep
91 : ! Number of MD steps held locally
92 :
93 : integer :: nshell3rd
94 : ! Number of shells for 3rd order IFC
95 :
96 : integer :: nshell4th
97 : ! Number of shells for 4th order IFC
98 :
99 : double precision :: U0
100 : ! Constant term contribution to energy
101 :
102 : double precision, allocatable :: Forces(:)
103 : ! Forces(3*natom*my_nstep)
104 : ! Forces computed from the fitted IFC and the displacements
105 :
106 : double precision, allocatable :: Phi1(:)
107 : ! Phi1(3*natom)
108 : ! First-order IFC
109 :
110 : type(Phi2_type) :: Phi2
111 : ! Second-order IFC
112 :
113 : double precision, allocatable :: Phi3(:,:,:,:)
114 : ! Phi3(3,3,3,Shell3at%nshell)
115 : ! Third-order IFC
116 :
117 : double precision, allocatable :: Phi4(:,:,:,:,:)
118 : ! Phi3(3,3,3,3,Shell4at%nshell)
119 : ! Fourth-order IFC
120 :
121 : double precision, allocatable :: Phi1Ui(:)
122 : ! Phi1Ui(my_nstep)
123 : ! First-order IFC contribution to energy
124 :
125 : double precision, allocatable :: Phi2UiUj(:)
126 : ! Phi2UiUj(my_nstep)
127 : ! Second-order IFC contribution to energy
128 :
129 : double precision, allocatable :: Phi3UiUjUk(:)
130 : ! Phi3UiUjUk(my_nstep)
131 : ! Third-order IFC contribution to energy
132 :
133 : double precision, allocatable :: Phi4UiUjUkUl(:)
134 : ! Phi4UiUjUkUl(my_nstep)
135 : ! Fourth-order IFC contribution to energy
136 :
137 : end type tdep_Model_type
138 :
139 : public :: tdep_init_phi2
140 : public :: tdep_destroy_phi2
141 : public :: tdep_model_init
142 : public :: tdep_model_free
143 : public :: tdep_calc_model
144 :
145 : contains
146 :
147 : !=====================================================================================================
148 :
149 46 : subroutine tdep_init_phi2(Phi2,dipdip,natom)
150 :
151 : type(Phi2_type),intent(out) :: Phi2
152 : logical, intent(in) :: dipdip
153 : integer, intent(in) :: natom
154 :
155 13524742 : ABI_CALLOC(Phi2%SR ,(3*natom,3*natom))
156 46 : if (dipdip) then
157 923727 : ABI_CALLOC(Phi2%Tot,(3*natom,3*natom))
158 923727 : ABI_CALLOC(Phi2%LR ,(3*natom,3*natom))
159 : end if
160 :
161 46 : end subroutine tdep_init_phi2
162 :
163 : !=====================================================================================================
164 :
165 46 : subroutine tdep_destroy_phi2(Phi2)
166 :
167 : type(Phi2_type),intent(inout) :: Phi2
168 :
169 46 : ABI_FREE(Phi2%SR)
170 46 : ABI_SFREE(Phi2%Tot)
171 46 : ABI_SFREE(Phi2%LR)
172 :
173 46 : end subroutine tdep_destroy_phi2
174 :
175 : !=====================================================================================================
176 :
177 44 : subroutine tdep_model_init(Model, Invar, Shell3at, Shell4at)
178 :
179 : type(tdep_Model_type),intent(out) :: Model
180 : type(atdep_dataset_type),intent(in) :: Invar
181 : type(Shell_type),intent(in) :: Shell3at, Shell4at
182 :
183 44 : Model%order = Invar%order
184 44 : Model%natom = Invar%natom
185 44 : Model%natom_unitcell = Invar%natom_unitcell
186 44 : Model%my_nstep = Invar%my_nstep
187 44 : Model%nshell3rd = 1
188 44 : Model%nshell4th = 1
189 :
190 44 : Model%U0 = zero
191 :
192 417588 : ABI_CALLOC(Model%Forces, (3*Model%natom*Model%my_nstep))
193 1084 : ABI_CALLOC(Model%Phi1Ui ,(Model%my_nstep))
194 1040 : ABI_CALLOC(Model%Phi2UiUj ,(Model%my_nstep))
195 1040 : ABI_CALLOC(Model%Phi3UiUjUk ,(Model%my_nstep))
196 1040 : ABI_CALLOC(Model%Phi4UiUjUkUl,(Model%my_nstep))
197 :
198 20778 : ABI_CALLOC(Model%Phi1,(3*Model%natom))
199 :
200 44 : call tdep_init_phi2(Model%Phi2,Invar%loto,Invar%natom)
201 :
202 44 : if (Model%order.ge.3) then
203 14 : Model%nshell3rd = Shell3at%nshell
204 2002 : ABI_CALLOC(Model%Phi3,(3,3,3,Model%nshell3rd))
205 : end if
206 :
207 44 : if (Model%order.eq.4) then
208 8 : Model%nshell4th = Shell4at%nshell
209 5590 : ABI_CALLOC(Model%Phi4,(3,3,3,3,Model%nshell4th))
210 : end if
211 :
212 44 : end subroutine tdep_model_init
213 :
214 : !=====================================================================================================
215 :
216 44 : subroutine tdep_model_free(Model)
217 :
218 : type(tdep_Model_type),intent(inout) :: Model
219 :
220 44 : ABI_FREE(Model%Forces)
221 44 : ABI_FREE(Model%Phi1Ui)
222 44 : ABI_FREE(Model%Phi2UiUj)
223 44 : ABI_FREE(Model%Phi3UiUjUk)
224 44 : ABI_FREE(Model%Phi4UiUjUkUl)
225 :
226 44 : ABI_FREE(Model%Phi1)
227 :
228 44 : call tdep_destroy_phi2(Model%Phi2)
229 :
230 44 : if (Model%order.ge.3) then
231 14 : ABI_FREE(Model%Phi3)
232 : end if
233 :
234 44 : if (Model%order.eq.4) then
235 8 : ABI_FREE(Model%Phi4)
236 : end if
237 :
238 44 : end subroutine tdep_model_free
239 :
240 : !====================================================================================================
241 :
242 44 : subroutine tdep_calc_model(Model,MD,Invar,MPIdata)
243 :
244 : type(tdep_Model_type),intent(inout) :: Model
245 : type(tdep_Sampling_type),intent(in) :: MD
246 : type(atdep_dataset_type),intent(in) :: Invar
247 : type(MPI_enreg_type), intent(in) :: MPIdata
248 :
249 : integer :: ii,jj,istep,iatom
250 : double precision :: Delta_F2,Delta_U,Delta_U2
251 : double precision :: sigma,U_1,U_2,U_3,U_4,UMD
252 44 : double precision, allocatable :: tmp(:),Phi_tot(:)
253 44 : double precision, allocatable :: U_MD(:),U_TDEP(:),weights_tot(:)
254 : integer :: ierr
255 :
256 : ! Compute the different contributions to total energy from the model
257 :
258 44 : write(Invar%stdout,*)' '
259 44 : write(Invar%stdout,*) '#############################################################################'
260 44 : write(Invar%stdout,*) '######################### Energies, errors,... #############################'
261 44 : write(Invar%stdout,*) '#############################################################################'
262 :
263 : ! Compute U0, U_TDEP, Delta_U and write them in the output file
264 44 : write(Invar%stdout,'(a)') ' Thermodynamic quantities and convergence parameters of THE MODEL,'
265 44 : write(Invar%stdout,'(a)') ' as a function of the step number (energies in eV/atom and forces in Ha/bohr) :'
266 44 : if (Invar%order.eq.4) then
267 8 : write(Invar%stdout,'(a)') ' <U_TDEP> = U_0 + U_1 + U_2 + U_3 + U_4'
268 8 : write(Invar%stdout,'(2a)') ' with U_0 = < U_MD - sum_i Phi1 ui - 1/2 sum_ij Phi2 ui uj ',&
269 16 : & '- 1/6 sum_ijk Phi3 ui uj uk - 1/24 sum_ijkl Phi4 ui uj uk ul >'
270 8 : write(Invar%stdout,'(a)') ' and U_1 = < sum_i Phi1 ui >'
271 8 : write(Invar%stdout,'(a)') ' and U_2 = < 1/2 sum_ij Phi2 ui uj >'
272 8 : write(Invar%stdout,'(a)') ' and U_3 = < 1/6 sum_ijk Phi3 ui uj uk >'
273 8 : write(Invar%stdout,'(a)') ' and U_4 = < 1/24 sum_ijkl Phi4 ui uj uk ul >'
274 36 : else if (Invar%order.eq.3) then
275 6 : write(Invar%stdout,'(a)') ' <U_TDEP> = U_0 + U_1 + U_2 + U_3'
276 6 : write(Invar%stdout,'(a)') ' with U_0 = < U_MD - sum_i Phi1 ui - 1/2 sum_ij Phi2 ui uj - 1/6 sum_ijk Phi3 ui uj uk >'
277 6 : write(Invar%stdout,'(a)') ' and U_1 = < sum_i Phi1 ui >'
278 6 : write(Invar%stdout,'(a)') ' and U_2 = < 1/2 sum_ij Phi2 ui uj >'
279 6 : write(Invar%stdout,'(a)') ' and U_3 = < 1/6 sum_ijk Phi3 ui uj uk >'
280 : else
281 30 : write(Invar%stdout,'(a)') ' <U_TDEP> = U_0 + U_1 + U_2'
282 30 : write(Invar%stdout,'(a)') ' with U_0 = < U_MD - sum_i Phi1 ui - 1/2 sum_ij Phi2 ui uj >'
283 30 : write(Invar%stdout,'(a)') ' and U_1 = < sum_i Phi1 ui >'
284 30 : write(Invar%stdout,'(a)') ' and U_2 = < 1/2 sum_ij Phi2 ui uj >'
285 : end if
286 44 : write(Invar%stdout,'(a)') ' Delta_U = < U_MD - U_TDEP > '
287 44 : write(Invar%stdout,'(a)') ' Delta_U2= (< (U_MD - U_TDEP)^2 >)**0.5 '
288 44 : write(Invar%stdout,'(a)') ' Delta_F2= (< (F_MD - F_TDEP)^2 >)**0.5 '
289 44 : write(Invar%stdout,'(a)') ' Sigma = (< (F_MD - F_TDEP)^2 >/<F_MD**2>)**0.5 '
290 44 : if (Invar%order.eq.4) then
291 8 : write(Invar%stdout,'(2a)') ' <U_MD> U_0 U_1 U_2 ',&
292 16 : & ' U_3 U_4 Delta_U Delta_U2 Delta_F2 Sigma'
293 36 : else if (Invar%order.eq.3) then
294 6 : write(Invar%stdout,'(2a)') ' <U_MD> U_0 U_1 U_2 ',&
295 12 : & ' U_3 Delta_U Delta_U2 Delta_F2 Sigma'
296 : else
297 30 : write(Invar%stdout,'(2a)') ' <U_MD> U_0 U_1 U_2 ',&
298 60 : & ' Delta_U Delta_U2 Delta_F2 Sigma'
299 : end if
300 :
301 : ! Compute eucledian distance for forces
302 528 : ABI_MALLOC(tmp,(11)) ; tmp(:) =0.d0
303 996 : do istep=1,MD%my_nstep
304 140148 : do iatom=1,MD%natom
305 557560 : do ii=1,3
306 417456 : jj = ii + 3*(iatom-1) + 3*MD%natom*(istep-1)
307 417456 : tmp(4)=tmp(4)+(MD%Forces(jj)-Model%Forces(jj))**2*MD%weights(istep)
308 556608 : tmp(5)=tmp(5)+MD%Forces(jj)**2*MD%weights(istep)
309 : end do
310 : end do
311 : end do
312 : ! Compute energies
313 1084 : ABI_MALLOC(U_TDEP, (MD%nstep_tot)) ; U_TDEP(:)=0.d0
314 1040 : ABI_MALLOC(U_MD, (MD%nstep_tot)) ; U_MD(:) =0.d0
315 1040 : ABI_MALLOC(weights_tot,(MD%nstep_tot)) ; weights_tot(:)=0.d0
316 1084 : ABI_MALLOC(Phi_tot, (MPIdata%my_nstep)); Phi_tot(:)=0.d0
317 996 : do istep=1,MD%my_nstep
318 952 : tmp(7) =tmp(7) +MD%etot(istep)*MD%weights(istep)
319 952 : tmp(10)=tmp(10)+Model%Phi1Ui(istep)*MD%weights(istep)
320 952 : tmp(6) =tmp(6) +Model%Phi2UiUj(istep)*MD%weights(istep)
321 952 : tmp(8) =tmp(8) +Model%Phi3UiUjUk(istep)*MD%weights(istep)
322 996 : tmp(11)=tmp(11)+Model%Phi4UiUjUkUl(istep)*MD%weights(istep)
323 : end do
324 44 : call xmpi_sum(tmp,MPIdata%comm_step,ierr)
325 44 : tmp(1) = tmp(7)-tmp(10)-tmp(6)-tmp(8)-tmp(11)
326 996 : Phi_tot(:)=tmp(1)+Model%Phi1Ui(:)+Model%Phi2UiUj(:)+Model%Phi3UiUjUk(:)+Model%Phi4UiUjUkUl(:)
327 : call xmpi_gatherv(Phi_tot,MD%my_nstep,U_TDEP,MPIdata%nstep_all,MPIdata%shft_step,&
328 44 : & MPIdata%master,MPIdata%comm_step,ierr)
329 : call xmpi_gatherv(MD%etot,MD%my_nstep,U_MD,MPIdata%nstep_all,MPIdata%shft_step,&
330 44 : & MPIdata%master,MPIdata%comm_step,ierr)
331 : call xmpi_gatherv(MD%weights,MD%my_nstep,weights_tot,MPIdata%nstep_all,MPIdata%shft_step,&
332 44 : & MPIdata%master,MPIdata%comm_step,ierr)
333 996 : do istep=1,MD%nstep_tot
334 952 : tmp(2) =tmp(2) + (U_MD(istep)-U_TDEP(istep)) * weights_tot(istep)
335 996 : tmp(9) =tmp(9) + (U_MD(istep)-U_TDEP(istep))**2 * weights_tot(istep)
336 : end do
337 44 : Model%U0 =tmp(1) /real(MD%natom)
338 44 : UMD =tmp(7) /real(MD%natom)
339 44 : U_1 =tmp(10)/real(MD%natom)
340 44 : U_2 =tmp(6) /real(MD%natom)
341 44 : U_3 =tmp(8) /real(MD%natom)
342 44 : U_4 =tmp(11)/real(MD%natom)
343 44 : Delta_U =tmp(2) /real(MD%natom)
344 44 : Delta_U2 =tmp(9) /real(MD%natom)
345 44 : Delta_F2 =tmp(4) /real(MD%natom*3)
346 44 : if (tmp(5).eq.0.d0) then
347 1 : sigma = 0.d0
348 : else
349 43 : sigma = dsqrt(tmp(4)/tmp(5))
350 : end if
351 44 : if (Invar%order.eq.4) then
352 8 : write(Invar%stdout,'(10(f12.5,5x))') UMD*Ha_eV,Model%U0*Ha_eV,U_1*Ha_eV,U_2*Ha_eV,U_3*Ha_eV,U_4*Ha_eV,&
353 16 : & Delta_U*Ha_eV,Delta_U2**0.5*Ha_eV,Delta_F2**0.5,sigma
354 36 : else if (Invar%order.eq.3) then
355 6 : write(Invar%stdout,'(9(f12.5,5x))') UMD*Ha_eV,Model%U0*Ha_eV,U_1*Ha_eV,U_2*Ha_eV,U_3*Ha_eV,&
356 12 : & Delta_U*Ha_eV,Delta_U2**0.5*Ha_eV,Delta_F2**0.5,sigma
357 : else
358 30 : write(Invar%stdout,'(8(f12.5,5x))') UMD*Ha_eV,Model%U0*Ha_eV,U_1*Ha_eV,U_2*Ha_eV,&
359 60 : & Delta_U*Ha_eV,Delta_U2**0.5*Ha_eV,Delta_F2**0.5,sigma
360 : endif
361 44 : ABI_FREE(tmp)
362 44 : write(Invar%stdout,'(a,1x,f12.5)') ' NOTE : in the harmonic and classical limit (T>>T_Debye), U_2=3/2*kB*T=',&
363 88 : & 3.d0/2.d0*kb_HaK*Ha_eV*Invar%temperature
364 :
365 : ! Write : i) (U_TDEP vs U_MD) in etotMDvsTDEP.dat
366 : ! ii) (Model%Forces vs MD%Forces) in fcartMDvsTDEP.dat
367 44 : write(Invar%stdout,'(a)') ' '
368 44 : write(Invar%stdout,'(a)') ' See the etotMDvsTDEP.dat & fcartMDvsTDEP.dat files'
369 44 : if (MPIdata%iam_master) then
370 44 : open(unit=32,file=trim(Invar%output_prefix)//'_etotMDvsTDEP.dat')
371 44 : open(unit=33,file=trim(Invar%output_prefix)//'_fcartMDvsTDEP.dat')
372 44 : write(32,'(a)') '# Istep U_MD(Ha) U_TDEP(Ha)'
373 44 : write(33,'(a)') '# Forces_MD(Ha/bohr) Forces_TDEP(Ha/bohr)'
374 996 : do istep=1,MD%nstep_tot
375 996 : write(32,'(i6,1x,2(f17.6,1x))') istep,U_MD(istep),U_TDEP(istep)
376 : end do
377 996 : do istep=1,MD%my_nstep
378 140148 : do iatom=1,MD%natom
379 557560 : do ii=1,3
380 417456 : write(33,'(2(f17.10,1x))') MD%Forces (ii+3*(iatom-1)+3*MD%natom*(istep-1)),&
381 974064 : & Model%Forces(ii+3*(iatom-1)+3*MD%natom*(istep-1))
382 : end do
383 : end do
384 : end do
385 44 : close(32)
386 44 : close(33)
387 : end if
388 44 : ABI_FREE(U_MD)
389 44 : ABI_FREE(U_TDEP)
390 44 : ABI_FREE(Phi_tot)
391 44 : ABI_FREE(weights_tot)
392 :
393 88 : end subroutine tdep_calc_model
394 :
395 : !====================================================================================================
396 :
397 0 : end module m_tdep_model
398 : !!***
|
