Line data Source code
1 : !!****m* ABINIT/m_predict_pimd
2 : !! NAME
3 : !! m_predict_pimd
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2010-2026 ABINIT group (GG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_predict_pimd
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_pimd
27 : use m_xmpi
28 : use m_results_img
29 :
30 : use defs_abitypes, only : MPI_type
31 : use m_geometry, only : mkradim, mkrdim
32 : use m_pimd_langevin, only : pimd_langevin_npt, pimd_langevin_nvt
33 : use m_pimd_nosehoover, only : pimd_nosehoover_npt, pimd_nosehoover_nvt
34 :
35 : implicit none
36 :
37 : private
38 : !!***
39 :
40 : public :: predict_pimd
41 : !!***
42 :
43 : contains
44 : !!***
45 :
46 : !!****f* ABINIT/predict_pimd
47 : !! NAME
48 : !! predict_pimd
49 : !!
50 : !! FUNCTION
51 : !! Predicts new positions in Path Integral Molecular Dynamics
52 : !! Given the positions at time t and t-dtion, an estimation of the velocities at time t,
53 : !! the forces and an estimation of the stress at time t, and an estimation of the cell at time t,
54 : !! computes in the Path Integral Molecular Dynamics framework the new positions at time t+dtion,
55 : !! computes self-consistently the velocities, the stress and the cell at time t and produces
56 : !! an estimation of the velocities, stress and new cell at time t+dtion
57 : !! No change of acell and rprim at present.
58 : !!
59 : !! INPUTS
60 : !! imgmov=gives the algorithm to be used for prediction of new set of images
61 : !! itimimage=time index for image propagation (itimimage+1 is to be predicted here)
62 : !! itimimage_eff=time index in the history
63 : !! mpi_enreg=MPI-parallelisation information
64 : !! natom=dimension of vel_timimage and xred_timimage
65 : !! nimage=number of images (treated by current proc)
66 : !! nimage_tot=total number of images
67 : !! ntimimage_stored=number of time steps stored in the history
68 : !! results_gs_timimage(ntimimage,nimage)=datastructure that hold all the history of previous computations.
69 : !! pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
70 : !! prtvolimg=printing volume
71 : !!
72 : !! OUTPUT
73 : !!
74 : !! SIDE EFFECTS
75 : !! results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
76 : !! results_img(:,:)%acell(3)
77 : !! at input, history of the values of acell for all images
78 : !! at output, the predicted values of acell for all images
79 : !! results_img(:,:)%results_gs
80 : !! at input, history of the values of energies and forces for all images
81 : !! results_img(:,:)%rprim(3,3)
82 : !! at input, history of the values of rprim for all images
83 : !! at output, the predicted values of rprim for all images
84 : !! results_img(:,:)%vel(3,natom)
85 : !! at input, history of the values of vel for all images
86 : !! at output, the predicted values of vel for all images
87 : !! results_img(:,:)%vel_cell(3,3)
88 : !! at input, history of the values of vel_cell for all images
89 : !! at output, the predicted values of vel_cell for all images
90 : !! results_img(:,:)%xred(3,natom)
91 : !! at input, history of the values of xred for all images
92 : !! at output, the predicted values of xred for all images
93 : !!
94 : !! SOURCE
95 :
96 35 : subroutine predict_pimd(imgmov,itimimage,itimimage_eff,mpi_enreg,natom,nimage,nimage_tot,&
97 35 : & ntimimage_stored,pimd_param,prtvolimg,results_img)
98 :
99 : !Arguments ------------------------------------
100 : !scalars
101 : integer,intent(in) :: imgmov,itimimage,itimimage_eff,natom
102 : integer,intent(in) :: nimage,nimage_tot,ntimimage_stored,prtvolimg
103 : type(MPI_type),intent(in) :: mpi_enreg
104 : type(pimd_type),intent(inout) :: pimd_param
105 : !arrays
106 : type(results_img_type) :: results_img(nimage,ntimimage_stored)
107 :
108 : !Local variables-------------------------------
109 : !scalars
110 : integer :: ierr,ii,itime,itime_next,itime_prev
111 : real(dp) :: volume
112 : !arrays
113 : real(dp) :: rprimd(3,3),rprimd_next(3,3),rprimd_prev(3,3),vel_cell(3,3),vel_cell_next(3,3)
114 35 : real(dp),allocatable :: mpibuf(:),mpibuffer(:,:,:),mpibuffer_all(:,:,:)
115 35 : real(dp),allocatable :: etotal(:),forces(:,:,:),stressin(:,:,:)
116 35 : real(dp),allocatable :: vel(:,:,:),vel_next(:,:,:)
117 35 : real(dp),allocatable :: xred(:,:,:),xred_next(:,:,:),xred_prev(:,:,:)
118 :
119 : ! *************************************************************************
120 :
121 : !############# Parallelism stuff 1 #######################
122 :
123 : !Parallelism over image: only one process per image of the cell
124 35 : if (mpi_enreg%me_cell==0) then
125 :
126 35 : itime=itimimage_eff
127 35 : itime_prev=itime-1;if (itime_prev<1) itime_prev=ntimimage_stored
128 :
129 35 : if (mpi_enreg%paral_img==0.or.mpi_enreg%me_img==0) then
130 140 : ABI_MALLOC(xred,(3,natom,nimage_tot))
131 105 : ABI_MALLOC(xred_prev,(3,natom,nimage_tot))
132 105 : ABI_MALLOC(xred_next,(3,natom,nimage_tot))
133 105 : ABI_MALLOC(etotal,(nimage_tot))
134 105 : ABI_MALLOC(forces,(3,natom,nimage_tot))
135 105 : ABI_MALLOC(stressin,(3,3,nimage_tot))
136 105 : ABI_MALLOC(vel,(3,natom,nimage_tot))
137 105 : ABI_MALLOC(vel_next,(3,natom,nimage_tot))
138 : end if
139 :
140 : ! Parallelism: Gather positions/forces/velocities/stresses/energy from all images
141 35 : if (mpi_enreg%paral_img==1) then
142 0 : ABI_MALLOC(mpibuffer,(12,natom+1,nimage))
143 0 : do ii=1,nimage
144 0 : mpibuffer(1:3 ,1:natom,ii)=results_img(ii,itime)%xred(1:3,1:natom)
145 0 : mpibuffer(4:6 ,1:natom,ii)=results_img(ii,itime_prev)%xred(1:3,1:natom)
146 0 : mpibuffer(7:9 ,1:natom,ii)=results_img(ii,itime)%results_gs%fcart(1:3,1:natom)
147 0 : mpibuffer(10:12,1:natom,ii)=results_img(ii,itime)%vel(1:3,1:natom)
148 0 : mpibuffer(1:6 ,natom+1,ii)=results_img(ii,itime)%results_gs%strten(1:6)
149 0 : mpibuffer(7:12 ,natom+1,ii)=zero
150 : end do
151 0 : if (mpi_enreg%me_img==0) then
152 0 : ABI_MALLOC(mpibuffer_all,(12,natom+1,nimage_tot))
153 : end if
154 0 : call gather_array_img(mpibuffer,mpibuffer_all,mpi_enreg,only_one_per_img=.true.,allgather=.false.)
155 0 : ABI_FREE(mpibuffer)
156 0 : if (mpi_enreg%me_img==0) then
157 0 : do ii=1,nimage_tot
158 0 : xred (1:3,1:natom,ii)=mpibuffer_all(1:3 ,1:natom,ii)
159 0 : xred_prev(1:3,1:natom,ii)=mpibuffer_all(4:6 ,1:natom,ii)
160 0 : forces (1:3,1:natom,ii)=mpibuffer_all(7:9 ,1:natom,ii)
161 0 : vel (1:3,1:natom,ii)=mpibuffer_all(10:12,1:natom,ii)
162 0 : stressin (1,1,ii) =mpibuffer_all(1,natom+1,ii)
163 0 : stressin (2,2,ii) =mpibuffer_all(2,natom+1,ii)
164 0 : stressin (3,3,ii) =mpibuffer_all(3,natom+1,ii)
165 0 : stressin (3,2,ii) =mpibuffer_all(4,natom+1,ii)
166 0 : stressin (3,1,ii) =mpibuffer_all(5,natom+1,ii)
167 0 : stressin (2,1,ii) =mpibuffer_all(6,natom+1,ii)
168 0 : stressin (2,3,ii)=stressin (3,2,ii)
169 0 : stressin (1,3,ii)=stressin (3,1,ii)
170 0 : stressin (1,2,ii)=stressin (2,1,ii)
171 : end do
172 0 : ABI_FREE(mpibuffer_all)
173 : end if
174 0 : ABI_MALLOC(mpibuf,(nimage))
175 0 : if (mpi_enreg%me_img/=0) then
176 0 : ABI_MALLOC(etotal,(0))
177 : end if
178 0 : do ii=1,nimage
179 0 : mpibuf(ii)=results_img(ii,itime)%results_gs%etotal
180 : end do
181 0 : call xmpi_gather(mpibuf,nimage,etotal,nimage,0,mpi_enreg%comm_img,ierr)
182 0 : ABI_FREE(mpibuf)
183 0 : if (mpi_enreg%me_img/=0) then
184 0 : ABI_FREE(etotal)
185 : end if
186 :
187 : ! No parallelism: simply copy positions/forces/velocities/stresses/energy
188 : else
189 240 : do ii=1,nimage
190 2245 : xred (:,:,ii)=results_img(ii,itime)%xred(:,:)
191 2245 : xred_prev(:,:,ii)=results_img(ii,itime_prev)%xred(:,:)
192 2245 : forces (:,:,ii)=results_img(ii,itime)%results_gs%fcart(:,:)
193 2245 : vel (:,:,ii)=results_img(ii,itime)%vel(:,:)
194 205 : etotal ( ii)=results_img(ii,itime)%results_gs%etotal
195 205 : stressin (1,1,ii)=results_img(ii,itime)%results_gs%strten(1)
196 205 : stressin (2,2,ii)=results_img(ii,itime)%results_gs%strten(2)
197 205 : stressin (3,3,ii)=results_img(ii,itime)%results_gs%strten(3)
198 205 : stressin (3,2,ii)=results_img(ii,itime)%results_gs%strten(4)
199 205 : stressin (3,1,ii)=results_img(ii,itime)%results_gs%strten(5)
200 205 : stressin (2,1,ii)=results_img(ii,itime)%results_gs%strten(6)
201 205 : stressin (2,3,ii)=stressin (3,2,ii)
202 205 : stressin (1,3,ii)=stressin (3,1,ii)
203 240 : stressin (1,2,ii)=stressin (2,1,ii)
204 : end do
205 : end if
206 :
207 : ! Parallelism over image: only one process does the job
208 35 : if (mpi_enreg%paral_img==0.or.mpi_enreg%me_img==0) then
209 :
210 : ! ############# PIMD MD algorithm #########################
211 :
212 : ! Some useful quantities about the cells (common to all images)
213 : ! Take acell and rprim from 1st image
214 35 : call mkrdim(results_img(1,itime)%acell,results_img(1,itime)%rprim,rprimd)
215 35 : call mkrdim(results_img(1,itime_prev)%acell,results_img(1,itime_prev)%rprim,rprimd_prev)
216 455 : vel_cell(:,:)=results_img(1,itime)%vel_cell(:,:)
217 :
218 : ! Compute the volume of the supercell
219 : volume=rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
220 : & rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
221 35 : & rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
222 35 : volume=abs(volume)
223 :
224 55 : select case(imgmov)
225 :
226 : case(9,10) !Langevin
227 :
228 55 : select case(pimd_param%optcell)
229 : case(0) !NVT
230 : call pimd_langevin_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
231 20 : & rprimd,stressin,nimage_tot,vel,vel_next,volume,xred,xred_next,xred_prev)
232 : case(2) !NPT
233 : call pimd_langevin_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
234 : & rprimd,rprimd_next,rprimd_prev,stressin,nimage_tot,vel,vel_next,vel_cell,&
235 20 : & vel_cell_next,volume,xred,xred_next,xred_prev)
236 : end select
237 :
238 : case(13) !Nose Hoover chains
239 :
240 35 : select case(pimd_param%optcell)
241 : case(0) !NVT
242 : call pimd_nosehoover_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
243 15 : & rprimd,stressin,nimage_tot,vel,vel_next,volume,xred,xred_next,xred_prev)
244 : case(2) !NPT
245 : call pimd_nosehoover_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
246 : & rprimd,rprimd_next,rprimd_prev,stressin,nimage_tot,vel,vel_next,vel_cell,&
247 15 : & vel_cell_next,volume,xred,xred_next,xred_prev)
248 : end select
249 :
250 : end select
251 :
252 : ! ############# Parallelism stuff 2 ########################
253 :
254 : end if ! mpi_enreg%me_img==0
255 :
256 : ! Parallelism: dispatch results
257 : ! The trick: use (9,natom) to store xred,xred_next,vel for all atoms
258 : ! use (9 ) to store rprimd_next
259 140 : ABI_MALLOC(mpibuffer,(12,natom+4,nimage))
260 35 : if (mpi_enreg%paral_img==1) then
261 0 : if (mpi_enreg%me_img==0) then
262 0 : ABI_MALLOC(mpibuffer_all,(12,natom+4,nimage_tot))
263 0 : do ii=1,nimage_tot
264 0 : mpibuffer_all(1:3,1:natom,ii)=xred_next(1:3,1:natom,ii)
265 0 : mpibuffer_all(4:6,1:natom,ii)=xred(1:3,1:natom,ii)
266 0 : mpibuffer_all(7:9,1:natom,ii)=vel_next(1:3,1:natom,ii)
267 0 : mpibuffer_all(10:12,1:natom,ii)=vel(1:3,1:natom,ii)
268 0 : mpibuffer_all(1:3,natom+1,ii)=rprimd_next(1:3,1)
269 0 : mpibuffer_all(4:6,natom+1,ii)=rprimd_next(1:3,2)
270 0 : mpibuffer_all(7:9,natom+1,ii)=rprimd_next(1:3,3)
271 0 : mpibuffer_all(1:3,natom+3,ii)=rprimd(1:3,1)
272 0 : mpibuffer_all(4:6,natom+3,ii)=rprimd(1:3,2)
273 0 : mpibuffer_all(7:9,natom+3,ii)=rprimd(1:3,3)
274 0 : mpibuffer_all(1:3,natom+2,ii)=vel_cell_next(1:3,1)
275 0 : mpibuffer_all(4:6,natom+2,ii)=vel_cell_next(1:3,2)
276 0 : mpibuffer_all(7:9,natom+2,ii)=vel_cell_next(1:3,3)
277 0 : mpibuffer_all(1:3,natom+2,ii)=vel_cell(1:3,1)
278 0 : mpibuffer_all(4:6,natom+2,ii)=vel_cell(1:3,2)
279 0 : mpibuffer_all(7:9,natom+2,ii)=vel_cell(1:3,3)
280 : end do
281 : end if
282 0 : call scatter_array_img(mpibuffer,mpibuffer_all,mpi_enreg)
283 0 : if (mpi_enreg%me_img==0) then
284 0 : ABI_FREE(mpibuffer_all)
285 : end if
286 : else
287 240 : do ii=1,nimage
288 2245 : mpibuffer(1:3,1:natom,ii)=xred_next(1:3,1:natom,ii)
289 2245 : mpibuffer(4:6,1:natom,ii)=xred(1:3,1:natom,ii)
290 2245 : mpibuffer(7:9,1:natom,ii)=vel_next(1:3,1:natom,ii)
291 2245 : mpibuffer(10:12,1:natom,ii)=vel(1:3,1:natom,ii)
292 820 : mpibuffer(1:3,natom+1,ii)=rprimd_next(1:3,1)
293 820 : mpibuffer(4:6,natom+1,ii)=rprimd_next(1:3,2)
294 820 : mpibuffer(7:9,natom+1,ii)=rprimd_next(1:3,3)
295 820 : mpibuffer(1:3,natom+3,ii)=rprimd(1:3,1)
296 820 : mpibuffer(4:6,natom+3,ii)=rprimd(1:3,2)
297 820 : mpibuffer(7:9,natom+3,ii)=rprimd(1:3,3)
298 820 : mpibuffer(1:3,natom+2,ii)=vel_cell_next(1:3,1)
299 820 : mpibuffer(4:6,natom+2,ii)=vel_cell_next(1:3,2)
300 820 : mpibuffer(7:9,natom+2,ii)=vel_cell_next(1:3,3)
301 820 : mpibuffer(1:3,natom+2,ii)=vel_cell(1:3,1)
302 820 : mpibuffer(4:6,natom+2,ii)=vel_cell(1:3,2)
303 855 : mpibuffer(7:9,natom+2,ii)=vel_cell(1:3,3)
304 : end do
305 : end if
306 :
307 35 : if (mpi_enreg%paral_img==0.or.mpi_enreg%me_img==0) then
308 35 : ABI_FREE(xred)
309 35 : ABI_FREE(xred_prev)
310 35 : ABI_FREE(xred_next)
311 35 : ABI_FREE(etotal)
312 35 : ABI_FREE(forces)
313 35 : ABI_FREE(stressin)
314 35 : ABI_FREE(vel)
315 35 : ABI_FREE(vel_next)
316 : end if
317 :
318 : else
319 0 : ABI_MALLOC(mpibuffer,(12,natom+4,nimage))
320 :
321 : end if ! mpi_enreg%me_cell==0
322 :
323 : !Send results to all procs treating the same image
324 35 : call xmpi_bcast(mpibuffer,0,mpi_enreg%comm_cell,ierr)
325 :
326 : !Store results in final place
327 35 : itime=itimimage_eff
328 35 : itime_prev=itime-1;if (itime_prev<1) itime_prev=ntimimage_stored
329 35 : itime_next=itime+1;if (itime_next>ntimimage_stored) itime_next=1
330 240 : do ii=1,nimage
331 2245 : results_img(ii,itime_next)%xred(1:3,1:natom)=mpibuffer(1:3,1:natom,ii)
332 2245 : results_img(ii,itime)%xred(1:3,1:natom)=mpibuffer(4:6,1:natom,ii)
333 2245 : results_img(ii,itime_next)%vel(1:3,1:natom)=mpibuffer(7:9,1:natom,ii)
334 2280 : results_img(ii,itime)%vel(1:3,1:natom)=mpibuffer(10:12,1:natom,ii)
335 : end do
336 35 : if (pimd_param%optcell/=0) then
337 0 : do ii=1,nimage
338 0 : rprimd(1:3,1)=mpibuffer(1:3,natom+1,ii)
339 0 : rprimd(1:3,2)=mpibuffer(4:6,natom+1,ii)
340 0 : rprimd(1:3,3)=mpibuffer(7:9,natom+1,ii)
341 0 : call mkradim(results_img(ii,itime_next)%acell,results_img(ii,itime_next)%rprim,rprimd)
342 0 : rprimd_prev(1:3,1)=mpibuffer(1:3,natom+2,ii)
343 0 : rprimd_prev(1:3,2)=mpibuffer(4:6,natom+2,ii)
344 0 : rprimd_prev(1:3,3)=mpibuffer(7:9,natom+2,ii)
345 0 : call mkradim(results_img(ii,itime)%acell,results_img(ii,itime)%rprim,rprimd_prev)
346 0 : results_img(ii,itime_next)%vel_cell(1:3,1)=mpibuffer(1:3,natom+3,ii)
347 0 : results_img(ii,itime_next)%vel_cell(1:3,2)=mpibuffer(4:6,natom+3,ii)
348 0 : results_img(ii,itime_next)%vel_cell(1:3,3)=mpibuffer(7:9,natom+3,ii)
349 0 : results_img(ii,itime)%vel_cell(1:3,1)=mpibuffer(1:3,natom+4,ii)
350 0 : results_img(ii,itime)%vel_cell(1:3,2)=mpibuffer(4:6,natom+4,ii)
351 0 : results_img(ii,itime)%vel_cell(1:3,3)=mpibuffer(7:9,natom+4,ii)
352 : end do
353 : end if
354 35 : ABI_FREE(mpibuffer)
355 :
356 35 : end subroutine predict_pimd
357 : !!***
358 :
359 : end module m_predict_pimd
360 : !!***
|
