Line data Source code
1 : !!****m* ABINIT/m_pimd_langevin
2 : !! NAME
3 : !! m_pimd_langevin
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2011-2026 ABINIT group (GG,MT)
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_pimd_langevin
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_pimd
28 : use m_random_zbq
29 :
30 : use m_matrix, only : matr3inv
31 : use m_geometry, only : xcart2xred, xred2xcart
32 :
33 :
34 : implicit none
35 :
36 : private
37 : !!***
38 :
39 : public :: pimd_langevin_npt
40 : public :: pimd_langevin_nvt
41 : !!***
42 :
43 : contains
44 : !!***
45 :
46 : !!****f* ABINIT/pimd_langevin_npt
47 : !! NAME
48 : !! pimd_langevin_npt
49 : !!
50 : !! FUNCTION
51 : !! Predicts new positions in Path Integral Molecular Dynamics using Langevin thermostat in the NPT ensemble.
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 : !!
58 : !! INPUTS
59 : !! etotal(trotter)=electronic total energy for all images
60 : !! itimimage=number of the current time for image propagation (itimimage+1 is to be predicted here)
61 : !! natom=dimension of vel_timimage and xred_timimage
62 : !! pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
63 : !! prtvolimg=printing volume
64 : !! rprimd(3,3)=dimensionless unit cell vectors (common to all images) at time t (present time step)
65 : !! rprimd_prev(3,3)=dimensionless unit cell vectors (common to all images) at time t-dt (previous time step)
66 : !! stressin(3,3,trotter)=electronic stress tensor for each image
67 : !! trotter=Trotter number (total number of images)
68 : !! volume=volume of unit cell (common to all images)
69 : !! xred(3,natom,trotter)=reduced coordinates of atoms for all images at time t (present time step)
70 : !! xred_prev(3,natom,trotter)=reduced coordinates of atoms for all images at time t-dt (previous time step)
71 : !!
72 : !! OUTPUT
73 : !! rprimd_next(3,3)=dimensionless unit cell vectors (common to all images) at time t+dt (next time step)
74 : !! xred_next(3,natom,trotter)=reduced coordinates of atoms for all images at time t+dt (next time step)
75 : !!
76 : !! SIDE EFFECTS
77 : !! forces(3,natom,trotter)=forces over atoms for all images
78 : !! at input, electronic forces
79 : !! at output, electronic forces + quantum spring contribution
80 : !! vel(3,natom,trotter)/vel_next(3,natom,trotter)/=velocies of atoms for all images
81 : !! at input, vel(itimimage) = values of the estimated vel at t
82 : !! at output, vel(itimimage) = values of the exact vel at t
83 : !! vel_next(itimimage) = values of the estimated vel at t+dt
84 : !! vel_cell(3,3)/vel_cell_next(3,3)=time derivative of cell parameters
85 : !! at input, vel_cell(itimimage) = values of the estimated cell vel at t
86 : !! at output, vel_cell(itimimage) = values of the exact cell vel at t
87 : !! vel_cell_next(itimimage) = values of the estimated cell vel at t+dt
88 : !!
89 : !! NOTES
90 : !! Here follows PIMD in the NPT ensemble within the Langevin barostat algorithm
91 : !! of Quigley and Probert: J. Chem. Phys. 120, 11432 (2004) [[cite:Quigley2004]]
92 : !! and Comput. Phys. Comm. 169, 322 (2005) [[cite:Quigley2005]]
93 : !!
94 : !! SOURCE
95 :
96 5 : subroutine pimd_langevin_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
97 5 : & rprimd,rprimd_next,rprimd_prev,stressin,trotter,vel,vel_next,vel_cell,&
98 5 : & vel_cell_next,volume,xred,xred_next,xred_prev)
99 :
100 : !Arguments ------------------------------------
101 : !scalars
102 : integer,intent(in) :: itimimage,natom,prtvolimg,trotter
103 : real(dp),intent(in) :: volume
104 : type(pimd_type),intent(in) :: pimd_param
105 : !arrays
106 : real(dp),intent(in) :: etotal(trotter),rprimd(3,3),rprimd_prev(3,3),stressin(3,3,trotter)
107 : real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
108 : real(dp),intent(out) :: rprimd_next(3,3),xred_next(3,natom,trotter)
109 : real(dp),intent(out) :: vel_next(3,natom,trotter),vel_cell_next(3,3)
110 : real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter),vel_cell(3,3)
111 :
112 : !Local variables-------------------------------
113 : !Options
114 : ! Option for the Langevin algorithm correction
115 : integer,parameter :: ilangevin=0
116 : ! The following option forces the total of forces to be zero
117 : ! It prevents the translation of the center of mass
118 : ! If it is zero, no constraint on mass center is applied
119 : integer,parameter :: zeroforce=1
120 : ! Tolerance for the SC cycle
121 : real(dp),parameter :: tolerance=tol9
122 :
123 : !scalars
124 : integer :: idum=-5
125 : integer :: constraint,iatom,ii,iimage,irestart,jj,ndof,prtstress
126 : real(dp) :: dtion,eharm,eharm2,epot,friction,frictionbar,initemp,kt,rescale_temp,scalebar
127 : real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol,tracepg,wg
128 : !arrays
129 : real, parameter :: identity(3,3)=reshape((/(one,(zero,ii=1,3),jj=1,2),one/),(/3,3/))
130 : real(dp) :: aleabar(3,3),constraint_output(2),ddh(3,3),diffstress(3,3)
131 : real(dp) :: dstrhh(3,3),fg(3,3),invrprimd(3,3)
132 : real(dp) :: langev_bar(3,3),pg(3,3),pgdh(3,3),stress_pimd(3,3,3),strtarget(6),tmp(3,3)
133 5 : real(dp),allocatable :: alea(:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:),forces_pimd_red(:,:)
134 5 : real(dp),allocatable :: fsup(:,:),hxredpoint(:,:,:),inertmass(:),langev(:,:)
135 5 : real(dp),allocatable :: mass(:,:),quantummass(:),spring(:,:)
136 5 : real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:)
137 5 : real(dp),allocatable :: xredpoint(:,:,:)
138 :
139 : ! *************************************************************************
140 :
141 : !############# Initializations ###########################
142 :
143 : !Allocation of local arrays
144 20 : ABI_MALLOC(xcart,(3,natom,trotter))
145 15 : ABI_MALLOC(xcart_prev,(3,natom,trotter))
146 15 : ABI_MALLOC(xcart_next,(3,natom,trotter))
147 15 : ABI_MALLOC(forces_orig,(3,natom,trotter))
148 15 : ABI_MALLOC(forces_pimd,(3,natom,trotter))
149 15 : ABI_MALLOC(inertmass,(natom))
150 10 : ABI_MALLOC(quantummass,(natom))
151 :
152 : !Fill in the local variables
153 5 : ndof=3*natom*trotter
154 5 : rescale_temp=one;if(zeroforce==1)rescale_temp=dble(ndof)/dble(ndof-3)
155 65 : quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass
156 65 : inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
157 5 : initemp=pimd_param%mdtemp(1)/rescale_temp
158 5 : thermtemp=pimd_param%mdtemp(2)
159 5 : dtion=pimd_param%dtion
160 5 : kt=thermtemp*kb_HaK
161 5 : friction=pimd_param%friction
162 5 : wg=pimd_param%bmass
163 35 : strtarget(:)=pimd_param%strtarget(:) ! imposed stress tensor
164 5 : frictionbar=pimd_param%frictionbar ! friction coeff of barostat
165 5 : scalebar=sqrt(two*frictionbar*wg*kt/dtion)
166 135 : forces_orig=forces
167 5 : constraint=0
168 :
169 : !Masses and spring constants
170 10 : ABI_MALLOC(mass,(natom,1))
171 10 : ABI_MALLOC(spring,(natom,1))
172 5 : call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,0,trotter)
173 :
174 : !Initialize random forces
175 15 : ABI_MALLOC(alea,(3,natom,trotter))
176 20 : ABI_MALLOC(langev,(natom,trotter))
177 35 : langev(:,1)=sqrt(two*friction*inertmass(:)*kt/dtion)
178 : if(ilangevin==1)then
179 : langev(:,1)=langev(:,1)*sqrt(one-(friction*dtion/(two*inertmass(:))))
180 : end if
181 :
182 : !Random number generator initialization
183 5 : if(itimimage<=1) then
184 1 : call pimd_langevin_random_init(pimd_param%irandom,idum)
185 : end if
186 :
187 : !Compute cartesian coordinates
188 10 : do iimage=1,trotter
189 5 : call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage))
190 10 : call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
191 : end do
192 :
193 : !Determine if it is a restart or not
194 : !If this is a calculation from scratch,generate random distribution of velocities
195 5 : irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel,vel_cell)
196 :
197 : !Initialize derivatives
198 5 : if (mod(irestart,10)==0) then
199 1 : call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
200 : end if
201 : !vel_cell does not depend on Trotter...
202 5 : ddh(:,:)=vel_cell(:,:);if (irestart<10) ddh=zero
203 :
204 5 : if (itimimage<=1) then
205 :
206 : ! ========================= FIRST TIME STEP =======================================
207 :
208 3 : ABI_MALLOC(hxredpoint,(3,natom,trotter))
209 3 : ABI_MALLOC(xredpoint,(3,natom,trotter))
210 3 : ABI_MALLOC(forces_pimd_red,(3,natom))
211 2 : ABI_MALLOC(fsup,(3,natom))
212 :
213 2 : do iimage=1,trotter
214 182 : hxredpoint(:,:,iimage)=vel(:,:,iimage) - matmul(ddh(:,:),xred(:,:,iimage))
215 : end do
216 1 : call matr3inv(rprimd,invrprimd)
217 2 : do iimage=1,trotter
218 266 : xredpoint(:,:,iimage)=matmul(invrprimd(:,:),hxredpoint(:,:,iimage))
219 : end do
220 :
221 : ! Generate random numbers
222 1 : call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce)
223 :
224 : ! Compute PIMD and Langevin contributions to forces
225 1 : call pimd_forces(forces,natom,spring,0,trotter,xcart)
226 1 : call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint)
227 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
228 1 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
229 40 : tmp=matmul(invrprimd,ddh)
230 52 : pg=wg*matmul(ddh,invrprimd)
231 1 : tracepg=pg(1,1)+pg(2,2)+pg(3,3)
232 :
233 : ! Taylor algorithm
234 2 : do iimage=1,trotter
235 1 : call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red)
236 157 : fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage))
237 8 : do iatom=1,natom
238 : xred_next(:,iatom,iimage)=xred(:,iatom,iimage)+dtion*xredpoint(:,iatom,iimage) + &
239 : & half*( &
240 : & forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom)- &
241 : & (tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) &
242 25 : & )*dtion*dtion/inertmass(iatom)
243 : end do
244 : end do
245 :
246 : ! predict rprimd at time t+dt from taylor algorithm
247 1 : call pimd_langevin_random_bar(aleabar,pimd_param%irandom,idum)
248 52 : langev_bar=matmul(aleabar,rprimd)*scalebar
249 1 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart)
250 1 : diffstress=pimd_diff_stress(stress_pimd,strtarget)
251 40 : dstrhh=matmul(diffstress,rprimd)
252 40 : pgdh=matmul(pg,ddh)
253 1 : temperature1=pimd_temperature(mass,hxredpoint)*rescale_temp
254 13 : fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature1*kb_HaK*rprimd(:,:)
255 :
256 : rprimd_next(:,:)=rprimd(:,:) + dtion*ddh(:,:) + half*( &
257 : & fg(:,:)-wg*frictionbar*ddh(:,:)+langev_bar &
258 13 : & )*dtion*dtion/wg
259 :
260 : ! Recompute xcart_next
261 2 : do iimage=1,trotter
262 2 : call xred2xcart(natom,rprimd_next,xcart_next(:,:,iimage),xred_next(:,:,iimage))
263 : end do
264 :
265 : ! Compute stress tensor at t from virial theorem
266 1 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
267 :
268 : ! Translate from pressure to stress for print
269 40 : stress_pimd=-stress_pimd
270 :
271 : ! Compute temperature at current step
272 1 : temperature1=pimd_temperature(mass,vel)*rescale_temp
273 :
274 : ! Estimate the velocities at t+dt
275 2 : do iimage=1,trotter
276 8 : do iatom=1,natom
277 : xredpoint(:,iatom,iimage)=xredpoint(:,iatom,iimage)+(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom)- &
278 25 : (tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))))*dtion/inertmass(iatom)
279 : end do
280 : end do
281 13 : ddh(:,:)=ddh(:,:)+( fg(:,:)-wg*frictionbar*ddh(:,:)+langev_bar(:,:) )*dtion/wg
282 2 : do iimage=1,trotter
283 2 : call xred2xcart(natom,rprimd_next,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
284 : end do
285 2 : do iimage=1,trotter
286 8 : do iatom=1,natom
287 115 : vel_next(:,iatom,iimage)=hxredpoint(:,iatom,iimage)+matmul(ddh(:,:),xred(:,iatom,iimage))
288 : end do
289 : end do
290 :
291 : ! Compute new temperature
292 1 : temperature2=pimd_temperature(mass,vel)*rescale_temp
293 :
294 26 : vel=xredpoint !scaled velocities transmitted to step 2
295 :
296 1 : ABI_FREE(hxredpoint)
297 1 : ABI_FREE(xredpoint)
298 1 : ABI_FREE(forces_pimd_red)
299 1 : ABI_FREE(fsup)
300 :
301 : else
302 :
303 : ! ========================= OTHER TIME STEPS ======================================
304 :
305 : ! Additional allocations
306 12 : ABI_MALLOC(hxredpoint,(3,natom,trotter))
307 12 : ABI_MALLOC(xredpoint,(3,natom,trotter))
308 12 : ABI_MALLOC(forces_pimd_red,(3,natom))
309 8 : ABI_MALLOC(fsup,(3,natom))
310 :
311 4 : ddh=vel_cell(:,:)
312 :
313 8 : do iimage=1,trotter
314 728 : hxredpoint(:,:,iimage)=vel(:,:,iimage) - matmul(ddh(:,:),xred(:,:,iimage))
315 : end do
316 :
317 : ! first estimation of ddh, pg and its trace:
318 4 : call matr3inv(rprimd,invrprimd)
319 208 : pg=wg*matmul(ddh,invrprimd)
320 4 : tracepg=pg(1,1)+pg(2,2)+pg(3,3)
321 :
322 : ! Momenta hxredpoint = H ds/dt: estimation
323 8 : do iimage=1,trotter
324 632 : hxredpoint(:,:,iimage)=matmul(rprimd,vel(:,:,iimage))
325 : end do
326 :
327 : ! Compute temperature at t
328 4 : temperature1=pimd_temperature(mass,hxredpoint)*rescale_temp
329 :
330 : ! Generate random numbers
331 4 : call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce)
332 :
333 : ! Generate random numbers for the barostat
334 4 : call pimd_langevin_random_bar(aleabar,pimd_param%irandom,idum)
335 208 : langev_bar=matmul(aleabar,rprimd)*scalebar
336 :
337 : ! Compute PIMD and Langevin contributions to forces
338 4 : call pimd_forces(forces,natom,spring,0,trotter,xcart)
339 4 : call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint)
340 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
341 4 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
342 :
343 : ! Compute difference between instantaneous stress and imposed stress (barostat)
344 4 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart)
345 :
346 4 : diffstress=pimd_diff_stress(stress_pimd,strtarget)
347 :
348 : ! Compute "force" on supercell vectors
349 160 : dstrhh=matmul(diffstress,rprimd)
350 160 : pgdh=matmul(pg,ddh)
351 52 : fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature1*kb_HaK*rprimd(:,:)
352 :
353 : ! Evolve the supercell (fist estimation)
354 52 : rprimd_next=two*rprimd-rprimd_prev+(fg-wg*frictionbar*ddh+langev_bar)*dtion*dtion/wg
355 :
356 : ! Evolve atomic positions (first estimation)
357 160 : tmp=matmul(invrprimd,ddh)
358 8 : do iimage=1,trotter
359 4 : call xcart2xred(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
360 628 : fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage))
361 4 : call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red)
362 :
363 32 : do iatom=1,natom
364 : xred_next(:,iatom,iimage)= &
365 : & two*xred(:,iatom,iimage) - xred_prev(:,iatom,iimage) &
366 : & +(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom) &
367 : & -tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) &
368 100 : & *dtion*dtion/inertmass(iatom)
369 : end do
370 : end do
371 :
372 : ! Self-consistent loop
373 4 : temperature2=pimd_temperature(mass,xredpoint)*rescale_temp
374 4 : temp2_prev=temperature2; tol=one
375 :
376 25 : do while (tol>tolerance)
377 : ! Reestimate dH/dt at t
378 273 : ddh(:,:)=(rprimd_next(:,:)-rprimd_prev(:,:))/(two*dtion)
379 :
380 : ! Reestimate the scaled velocities at t
381 42 : do iimage=1,trotter
382 525 : xredpoint(:,:,iimage)=(xred_next(:,:,iimage)-xred_prev(:,:,iimage))/(two*dtion)
383 42 : call xred2xcart(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
384 : end do
385 : ! Reestimate the forces
386 21 : call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint)
387 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
388 21 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
389 : ! Compute variation of temperature (to check convergence of SC loop)
390 21 : temperature2=pimd_temperature(mass,xredpoint)*rescale_temp
391 21 : tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev)
392 21 : temp2_prev=temperature2
393 : ! Recompute the temperature
394 21 : temperature2=pimd_temperature(mass,hxredpoint)*rescale_temp
395 : ! Recompute pg
396 1092 : pg=wg*matmul(ddh,invrprimd)
397 21 : tracepg=pg(1,1)+pg(2,2)+pg(3,3)
398 : ! Recompute difference between instantaneous stress and imposed stress (barostat)
399 21 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart)
400 :
401 21 : diffstress=pimd_diff_stress(stress_pimd,strtarget)
402 :
403 : ! Recompute "force" on supercell vectors
404 840 : dstrhh=matmul(diffstress,rprimd)
405 840 : pgdh=matmul(pg,ddh)
406 273 : fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature2*kb_HaK*rprimd(:,:)
407 : ! Evolve the supercell (better estimation)
408 273 : rprimd_next=two*rprimd-rprimd_prev+(fg-wg*frictionbar*ddh+langev_bar)*dtion*dtion/wg
409 :
410 : ! Evolve atomic positions (better estimation):
411 840 : tmp=matmul(invrprimd,ddh)
412 46 : do iimage=1,trotter
413 21 : call xcart2xred(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage))
414 3297 : fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage))
415 21 : call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red)
416 168 : do iatom=1,natom
417 : xred_next(:,iatom,iimage)= &
418 : & two*xred(:,iatom,iimage) - xred_prev(:,iatom,iimage) &
419 : & +(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom) &
420 : & -(tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof)))) &
421 525 : & *dtion*dtion/inertmass(iatom)
422 : end do
423 : end do
424 : end do ! End self-consistent loop
425 :
426 : ! Computation of true temperature from true velocities at t
427 8 : do iimage=1,trotter
428 512 : vel(:,:,iimage)=hxredpoint(:,:,iimage)+matmul(ddh,xred(:,:,iimage))
429 : end do
430 4 : temperature2=pimd_temperature(mass,vel)*rescale_temp
431 :
432 : ! Computation of the real stress tensor at t
433 4 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
434 :
435 : ! translate from pressure to stress
436 160 : stress_pimd=-stress_pimd
437 :
438 : ! Deallocations (Verlet algo)
439 4 : ABI_FREE(xredpoint)
440 4 : ABI_FREE(hxredpoint)
441 4 : ABI_FREE(forces_pimd_red)
442 4 : ABI_FREE(fsup)
443 :
444 : end if ! itimimage==1
445 :
446 : !############# Final operations ############################
447 :
448 : !Compute contributions to energy
449 5 : call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring,trotter,xcart)
450 :
451 : !Print messages
452 5 : prtstress=1
453 : call pimd_print(constraint,constraint_output,&
454 : & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
455 : & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
456 : & stress_pimd,temperature2,&
457 5 : & pimd_param%traj_unit,trotter,vel,ddh,xcart,xred)
458 :
459 : !Compute cartesian coordinates
460 10 : do iimage=1,trotter
461 10 : call xred2xcart(natom,rprimd_next,xcart_next(:,:,iimage),xred_next(:,:,iimage))
462 : end do
463 :
464 : !If possible, estimate the velocities at t+dt
465 5 : if (itimimage>1) then
466 4 : call pimd_predict_vel(dtion,itimimage,natom,trotter,0,xcart,xcart_next,xcart_prev,vel_next)
467 4 : call pimd_predict_velcell(dtion,itimimage,rprimd,rprimd_next,rprimd_prev,ddh)
468 : end if
469 :
470 : !Return cell velocities (does not depend on Trotter)
471 5 : vel_cell_next(:,:)=ddh(:,:)
472 :
473 : !Free memory
474 5 : ABI_FREE(xcart)
475 5 : ABI_FREE(xcart_prev)
476 5 : ABI_FREE(xcart_next)
477 5 : ABI_FREE(forces_orig)
478 5 : ABI_FREE(forces_pimd)
479 5 : ABI_FREE(inertmass)
480 5 : ABI_FREE(quantummass)
481 5 : ABI_FREE(mass)
482 5 : ABI_FREE(spring)
483 5 : ABI_FREE(alea)
484 5 : ABI_FREE(langev)
485 :
486 10 : end subroutine pimd_langevin_npt
487 : !!***
488 :
489 : !!****f* ABINIT/pimd_langevin_nvt
490 : !! NAME
491 : !! pimd_langevin_nvt
492 : !!
493 : !! FUNCTION
494 : !! Predicts new positions in Path Integral Molecular Dynamics using Langevin thermostat in the NVT ensemble.
495 : !! Given the positions at time t and t-dtion, an estimation of the velocities at time t,
496 : !! the forces and an estimation of the stress at time t, and an estimation of the cell at time t,
497 : !! computes in the Path Integral Molecular Dynamics framework the new positions at time t+dtion,
498 : !! computes self-consistently the velocities, the stress and the cell at time t and produces
499 : !! an estimation of the velocities, stress and new cell at time t+dtion
500 : !!
501 : !! INPUTS
502 : !! etotal(trotter)=electronic total energy for all images
503 : !! itimimage=number of the current time for image propagation (itimimage+1 is to be predicted here)
504 : !! natom=dimension of vel_timimage and xred_timimage
505 : !! pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
506 : !! prtvolimg=printing volume
507 : !! rprimd(3,3)=dimensionless unit cell vectors (common to all images)
508 : !! stressin(3,3,trotter)=electronic stress tensor for each image
509 : !! trotter=Trotter number (total number of images)
510 : !! volume=volume of unit cell (common to all images)
511 : !! xred(3,natom,trotter)=reduced coordinates of atoms for all images at time t (present time step)
512 : !! xred_prev(3,natom,trotter)=reduced coordinates of atoms for all images at time t-dt (previous time step)
513 : !!
514 : !! OUTPUT
515 : !! xred_next(3,natom,trotter)=reduced coordinates of atoms for all images at time t+dt (next time step)
516 : !!
517 : !! SIDE EFFECTS
518 : !! forces(3,natom,trotter)=forces over atoms for all images
519 : !! at input, electronic forces
520 : !! at output, electronic forces + quantum spring contribution
521 : !! vel(3,natom,trotter)/vel_next(3,natom,trotter)/=velocies of atoms for all images
522 : !! at input, vel(itimimage) = values of the estimated vel at t
523 : !! at output, vel(itimimage) = values of the exact vel at t
524 : !! vel_next(itimimage) = values of the estimated vel at t+dt
525 : !!
526 : !! NOTES
527 : !! See Quigley,Probert, JCP 120, 11432 (2004) [[cite:Quigley2004]], part III
528 : !!
529 : !! SOURCE
530 :
531 25 : subroutine pimd_langevin_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,&
532 25 : & rprimd,stressin,trotter,vel,vel_next,volume,xred,xred_next,xred_prev)
533 :
534 : !Arguments ------------------------------------
535 : !scalars
536 : integer,intent(in) :: itimimage,natom,prtvolimg,trotter
537 : real(dp),intent(in) :: volume
538 : type(pimd_type),intent(in) :: pimd_param
539 : !arrays
540 : real(dp),intent(in) :: etotal(trotter),rprimd(3,3),stressin(3,3,trotter)
541 : real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter)
542 : real(dp),intent(out) :: xred_next(3,natom,trotter),vel_next(3,natom,trotter)
543 : real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter)
544 :
545 : !Local variables-------------------------------
546 : !Options
547 : ! Option for the Langevin algorithm correction
548 : integer,parameter :: ilangevin=0
549 : ! Tolerance for the SC cycle
550 : real(dp),parameter :: tolerance=tol9
551 :
552 : !scalars
553 : integer :: idum=-5
554 : integer :: iimage,irestart,ndof,pitransform,prtstress,use_qtb,zeroforce
555 : real(dp) :: dtion,eharm,eharm2,epot,friction,initemp,kt,kt_,rescale_temp
556 : real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol
557 : !arrays
558 50 : real(dp) :: constraint_output(2),spring_prim(natom),stress_pimd(3,3,3),vel_cell(3,3)
559 25 : real(dp),allocatable :: alea(:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:)
560 25 : real(dp),allocatable :: inertmass(:),langev(:,:),mass(:,:),quantummass(:),spring(:,:)
561 25 : real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:)
562 :
563 : ! *************************************************************************
564 :
565 25 : if (pimd_param%use_qtb==1.and.pimd_param%qtb_file_unit<=0) then
566 0 : ABI_BUG('piqtb_force not open!')
567 : end if
568 :
569 : !############# Initializations ###########################
570 :
571 25 : pitransform=pimd_param%pitransform
572 :
573 : !The following option forces the total of forces to be zero
574 : !It prevents the translation of the center of mass
575 : !If it is zero, no constraint on mass center is applied
576 25 : zeroforce=1
577 25 : if(pitransform==1) zeroforce=0
578 25 : if(pitransform==2) zeroforce=0
579 25 : if(pimd_param%constraint==1) zeroforce=0
580 :
581 : !Allocation of local arrays
582 100 : ABI_MALLOC(xcart,(3,natom,trotter))
583 75 : ABI_MALLOC(xcart_prev,(3,natom,trotter))
584 75 : ABI_MALLOC(xcart_next,(3,natom,trotter))
585 75 : ABI_MALLOC(forces_orig,(3,natom,trotter))
586 75 : ABI_MALLOC(forces_pimd,(3,natom,trotter))
587 75 : ABI_MALLOC(inertmass,(natom))
588 50 : ABI_MALLOC(quantummass,(natom))
589 :
590 : !Fill in the local variables
591 25 : use_qtb=pimd_param%use_qtb
592 25 : ndof=3*natom*trotter
593 25 : rescale_temp=one
594 25 : if(zeroforce==1) rescale_temp=dble(ndof)/dble(ndof-3)
595 205 : quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass
596 205 : inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass
597 25 : if(pitransform==1) inertmass=quantummass !compulsory for good definition of normal mode masses
598 55 : if(pitransform==2) inertmass=quantummass !compulsory for good definition of staging masses
599 25 : initemp=pimd_param%mdtemp(1)/rescale_temp
600 25 : thermtemp=pimd_param%mdtemp(2)
601 25 : friction=pimd_param%friction
602 25 : dtion=pimd_param%dtion
603 25 : kt=thermtemp*kb_HaK
604 1610 : forces_orig=forces
605 :
606 : !Masses and spring constants
607 15 : select case(pitransform)
608 : case(0)
609 30 : ABI_MALLOC(mass,(natom,1)) ! This second dimension is needed
610 30 : ABI_MALLOC(spring,(natom,1))
611 30 : ABI_MALLOC(langev,(natom,1))
612 : case(1,2)
613 40 : ABI_MALLOC(mass,(natom,trotter))
614 30 : ABI_MALLOC(spring,(natom,trotter))
615 55 : ABI_MALLOC(langev,(natom,trotter))
616 : end select
617 115 : spring_prim(:)=quantummass(:)*dble(trotter)*kt*kt
618 25 : call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,pitransform,trotter)
619 :
620 : !Initialize random forces
621 75 : ABI_MALLOC(alea,(3,natom,trotter))
622 25 : if (use_qtb==0) then
623 290 : langev(:,:)=sqrt(two*friction*mass(:,:)*kt/dtion)
624 : else
625 0 : langev(:,:)=sqrt(two*friction*mass(:,:))
626 : end if
627 :
628 : !Random number generator initialization
629 25 : if(itimimage<=1) then
630 4 : call pimd_langevin_random_init(pimd_param%irandom,idum)
631 : end if
632 :
633 : !Compute cartesian coordinates
634 145 : do iimage=1,trotter
635 120 : call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage))
636 145 : call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage))
637 : end do
638 :
639 : !Determine if it is a restart or not
640 : !If this is a calculation from scratch,generate random distribution of velocities
641 25 : irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel)
642 25 : if (irestart==0) then
643 4 : call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon)
644 : end if
645 :
646 : !Compute temperature at t
647 25 : temperature1=pimd_temperature(mass,vel)*rescale_temp
648 :
649 : !################## Images evolution #####################
650 :
651 : !Generate random numbers
652 25 : if (use_qtb==0) then
653 25 : call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce)
654 : else
655 0 : call pimd_langevin_random_qtb(alea,langev,mass,natom,pimd_param%qtb_file_unit,trotter,zeroforce)
656 : end if
657 :
658 : !Compute PIMD and Langevin contributions to forces
659 25 : call pimd_coord_transform(xcart,1,natom,pitransform,trotter)
660 25 : call pimd_force_transform(forces,1,natom,pitransform,trotter) !compute staging forces
661 25 : call pimd_forces(forces,natom,spring,pitransform,trotter,xcart)
662 25 : call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,vel)
663 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
664 25 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
665 :
666 : !Compute atomic positions at t+dt
667 25 : if (itimimage<=1) then
668 :
669 : ! === 1st time step: single Taylor algorithm
670 : ! Predict positions
671 : call pimd_predict_taylor(dtion,forces_pimd,mass,natom,trotter,&
672 4 : & vel,xcart,xcart_next)
673 :
674 : ! Compute new temperature
675 4 : temperature2=pimd_temperature(mass,vel)*rescale_temp
676 :
677 : else
678 :
679 : ! === Other time steps: Verlet algorithm + SC cycle
680 : ! Predict positions
681 21 : call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter)
682 : call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
683 21 : & xcart,xcart_next,xcart_prev)
684 : ! Self-consistent loop
685 21 : temperature2=pimd_temperature(mass,vel)*rescale_temp
686 21 : temp2_prev=temperature2; tol=one
687 106 : do while (tol>tolerance)
688 : ! Recompute a (better) estimation of the velocity at time step t
689 5702 : vel = (xcart_next - xcart_prev) / (two*dtion)
690 85 : temperature2=pimd_temperature(mass,vel)*rescale_temp
691 : ! Reestimate the force
692 : call pimd_langevin_forces(alea,forces,forces_pimd,friction,&
693 85 : & langev,mass,natom,trotter,vel)
694 : call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,&
695 85 : & mass,natom,trotter,pimd_param%wtatcon,xcart)
696 : ! Compute new positions
697 : call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,&
698 85 : & xcart,xcart_next,xcart_prev)
699 :
700 : ! Compute variation of temperature (to check convergence of SC loop)
701 85 : tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev)
702 85 : temp2_prev=temperature2
703 :
704 : end do ! End self-consistent loop
705 :
706 : end if ! itimimage==1
707 :
708 25 : call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter)
709 25 : call pimd_coord_transform(xcart,-1,natom,pitransform,trotter)
710 25 : call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter)
711 :
712 : !Compute contributions to energy
713 25 : call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring_prim,trotter,xcart)
714 :
715 : !Compute stress tensor at t
716 25 : if (use_qtb==0) then
717 25 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart)
718 : else
719 0 : call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,temperature2,thermtemp,trotter,vel,volume,xcart)
720 : end if
721 1000 : stress_pimd=-stress_pimd ! Translate pressure to stress
722 :
723 : !############# Final operations ############################
724 :
725 : !Print messages
726 25 : vel_cell=zero;prtstress=1;if (prtvolimg>=2) prtstress=0
727 25 : kt_=kt;if (use_qtb==1) kt_=temperature2*kb_HaK
728 : call pimd_print(pimd_param%constraint,constraint_output,&
729 : & eharm,eharm2,epot,forces_pimd,inertmass,irestart,&
730 : & itimimage,kt_,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,&
731 : & stress_pimd,temperature2,&
732 25 : & pimd_param%traj_unit,trotter,vel,vel_cell,xcart,xred)
733 :
734 : !If possible, estimate the (transformed) velocities at t+dt
735 25 : call pimd_predict_vel(dtion,itimimage,natom,trotter,pitransform,xcart,xcart_next,xcart_prev,vel_next)
736 :
737 : !Come back to reduced coordinates
738 145 : do iimage=1,trotter
739 145 : call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage))
740 : end do
741 :
742 : !Free memory
743 25 : ABI_FREE(xcart)
744 25 : ABI_FREE(xcart_prev)
745 25 : ABI_FREE(xcart_next)
746 25 : ABI_FREE(forces_orig)
747 25 : ABI_FREE(forces_pimd)
748 25 : ABI_FREE(inertmass)
749 25 : ABI_FREE(quantummass)
750 25 : ABI_FREE(mass)
751 25 : ABI_FREE(spring)
752 25 : ABI_FREE(alea)
753 25 : ABI_FREE(langev)
754 :
755 25 : end subroutine pimd_langevin_nvt
756 : !!***
757 :
758 : end module m_pimd_langevin
759 : !!***
|
