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