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