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