TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 2011-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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=voume 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,trotter)=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) and Comput. Phys. Comm. 169, 322 (2005)
PARENTS
predict_pimd
CHILDREN
matr3inv,pimd_apply_constraint,pimd_energies,pimd_forces,pimd_initvel pimd_langevin_forces,pimd_langevin_random,pimd_langevin_random_bar pimd_langevin_random_init,pimd_mass_spring,pimd_print,pimd_stresses xcart2xred,xred2xcart
SOURCE
66 #if defined HAVE_CONFIG_H 67 #include "config.h" 68 #endif 69 70 #include "abi_common.h" 71 72 subroutine pimd_langevin_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,& 73 & rprimd,rprimd_next,rprimd_prev,stressin,trotter,vel,vel_cell,& 74 & volume,xred,xred_next,xred_prev) 75 76 77 use defs_basis 78 use m_profiling_abi 79 use m_pimd 80 use m_random_zbq 81 82 use m_geometry, only : xcart2xred, xred2xcart 83 84 !This section has been created automatically by the script Abilint (TD). 85 !Do not modify the following lines by hand. 86 #undef ABI_FUNC 87 #define ABI_FUNC 'pimd_langevin_npt' 88 use interfaces_32_util 89 !End of the abilint section 90 91 implicit none 92 93 !Arguments ------------------------------------ 94 !scalars 95 integer,intent(in) :: itimimage,natom,prtvolimg,trotter 96 real(dp),intent(in) :: volume 97 type(pimd_type),intent(in) :: pimd_param 98 !arrays 99 real(dp),intent(in) :: etotal(trotter),rprimd(3,3),rprimd_prev(3,3),stressin(3,3,trotter) 100 real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter) 101 real(dp),intent(out) :: rprimd_next(3,3),xred_next(3,natom,trotter) 102 real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter),vel_cell(3,3,trotter) 103 104 !Local variables------------------------------- 105 !Options 106 ! Option for the Langevin algorithm correction 107 integer,parameter :: ilangevin=0 108 ! The following option forces the total of forces to be zero 109 ! It prevents the translation of the center of mass 110 ! If it is zero, no constraint on mass center is applied 111 integer,parameter :: zeroforce=1 112 ! Tolerance for the SC cycle 113 real(dp),parameter :: tolerance=tol7 114 115 !scalars 116 integer :: idum=-5 117 integer :: constraint,iatom,ii,iimage,irestart,jj,ndof,prtstress 118 real(dp) :: dtion,eharm,eharm2,epot,friction,frictionbar,initemp,kt,rescale_temp,scalebar 119 real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol,tracepg,wg 120 !arrays 121 real, parameter :: identity(3,3)=reshape((/(one,(zero,ii=1,3),jj=1,2),one/),(/3,3/)) 122 real(dp) :: aleabar(3,3),constraint_output(2),ddh(3,3),diffstress(3,3) 123 real(dp) :: dstrhh(3,3),fg(3,3),invrprimd(3,3) 124 real(dp) :: langev_bar(3,3),pg(3,3),pgdh(3,3),stress_pimd(3,3,3),strtarget(6),tmp(3,3) 125 real(dp),allocatable :: alea(:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:),forces_pimd_red(:,:) 126 real(dp),allocatable :: fsup(:,:),hxredpoint(:,:,:),inertmass(:),langev(:,:) 127 real(dp),allocatable :: mass(:,:),quantummass(:),spring(:,:) 128 real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:) 129 real(dp),allocatable :: xredpoint(:,:,:) 130 131 ! ************************************************************************* 132 133 !############# Initializations ########################### 134 135 !Allocation of local arrays 136 ABI_ALLOCATE(xcart,(3,natom,trotter)) 137 ABI_ALLOCATE(xcart_prev,(3,natom,trotter)) 138 ABI_ALLOCATE(xcart_next,(3,natom,trotter)) 139 ABI_ALLOCATE(forces_orig,(3,natom,trotter)) 140 ABI_ALLOCATE(forces_pimd,(3,natom,trotter)) 141 ABI_ALLOCATE(inertmass,(natom)) 142 ABI_ALLOCATE(quantummass,(natom)) 143 144 !Fill in the local variables 145 ndof=3*natom*trotter 146 rescale_temp=one;if(zeroforce==1)rescale_temp=dble(ndof)/dble(ndof-3) 147 quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass 148 inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass 149 initemp=pimd_param%mdtemp(1)/rescale_temp; thermtemp=pimd_param%mdtemp(2) 150 dtion=pimd_param%dtion 151 kt=thermtemp*kb_HaK 152 friction=pimd_param%vis 153 wg=pimd_param%bmass 154 strtarget(:)=pimd_param%strtarget(:) ! imposed stress tensor 155 frictionbar=pimd_param%friction ! friction coeff of barostat 156 scalebar=sqrt(two*frictionbar*wg*kt/dtion) 157 forces_orig=forces 158 159 !Masses and spring constants 160 ABI_ALLOCATE(mass,(natom,1)) 161 ABI_ALLOCATE(spring,(natom,1)) 162 call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,0,trotter) 163 164 !Initialize random forces 165 ABI_ALLOCATE(alea,(3,natom,trotter)) 166 ABI_ALLOCATE(langev,(natom,trotter)) 167 langev(:,1)=sqrt(two*friction*inertmass(:)*kt/dtion) 168 if(ilangevin==1)then 169 langev(:,1)=langev(:,1)*sqrt(one-(friction*dtion/(two*inertmass(:)))) 170 end if 171 172 !Random number generator initialization 173 if(itimimage<=1) then 174 call pimd_langevin_random_init(pimd_param%irandom,idum) 175 end if 176 177 !Compute cartesian coordinates 178 do iimage=1,trotter 179 call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage)) 180 call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage)) 181 end do 182 183 !Determine if it is a restart or not 184 !If this is a calculation from scratch,generate random distribution of velocities 185 irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel,vel_cell) 186 187 !Initialize derivatives 188 if (mod(irestart,10)==0) then 189 call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon) 190 end if 191 !vel_cell does not depend on Trotter... 192 ddh=vel_cell(:,:,1);if (irestart<10) ddh=zero 193 194 if (itimimage<=1) then 195 196 ! ========================= FIRST TIME STEP ======================================= 197 198 ABI_ALLOCATE(hxredpoint,(3,natom,trotter)) 199 ABI_ALLOCATE(xredpoint,(3,natom,trotter)) 200 ABI_ALLOCATE(forces_pimd_red,(3,natom)) 201 ABI_ALLOCATE(fsup,(3,natom)) 202 203 do iimage=1,trotter 204 hxredpoint(:,:,iimage)=vel(:,:,iimage) - matmul(ddh(:,:),xred(:,:,iimage)) 205 end do 206 call matr3inv(rprimd,invrprimd) 207 do iimage=1,trotter 208 xredpoint(:,:,iimage)=matmul(invrprimd(:,:),hxredpoint(:,:,iimage)) 209 end do 210 211 ! Generate random numbers 212 call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce) 213 214 ! Compute PIMD and Langevin contributions to forces 215 call pimd_forces(forces,natom,spring,0,trotter,xcart) 216 call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint) 217 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 218 & mass,natom,trotter,pimd_param%wtatcon,xcart) 219 tmp=matmul(invrprimd,ddh) 220 pg=wg*matmul(ddh,invrprimd) 221 tracepg=pg(1,1)+pg(2,2)+pg(3,3) 222 223 ! Taylor algorithm 224 do iimage=1,trotter 225 call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red) 226 fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage)) 227 do iatom=1,natom 228 xred_next(:,iatom,iimage)=xred(:,iatom,iimage)+dtion*xredpoint(:,iatom,iimage) + & 229 & half*( & 230 & forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom)- & 231 & (tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) & 232 & )*dtion*dtion/inertmass(iatom) 233 end do 234 end do 235 236 ! predict rprimd at time t+dt from taylor algorithm 237 call pimd_langevin_random_bar(aleabar,pimd_param%irandom,idum) 238 langev_bar=matmul(aleabar,rprimd)*scalebar 239 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart) 240 diffstress=pimd_diff_stress(stress_pimd,strtarget) 241 dstrhh=matmul(diffstress,rprimd) 242 pgdh=matmul(pg,ddh) 243 temperature1=pimd_temperature(mass,hxredpoint)*rescale_temp 244 fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature1*kb_HaK*rprimd(:,:) 245 246 rprimd_next(:,:)=rprimd(:,:) + dtion*ddh(:,:) + half*( & 247 & fg(:,:)-wg*frictionbar*ddh(:,:)+langev_bar & 248 & )*dtion*dtion/wg 249 250 ! Recompute xcart_next 251 do iimage=1,trotter 252 call xred2xcart(natom,rprimd_next,xcart_next(:,:,iimage),xred_next(:,:,iimage)) 253 end do 254 255 ! Compute stress tensor at t from virial theorem 256 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart) 257 258 ! Translate from pressure to stress for print 259 stress_pimd=-stress_pimd 260 261 ! Compute temperature at current step 262 temperature1=pimd_temperature(mass,vel)*rescale_temp 263 264 ! Estimate the velocities at t+dt/2 265 vel=(xcart_next-xcart)/dtion 266 ddh=(rprimd_next-rprimd)/dtion 267 268 ! Compute new temperature 269 temperature2=pimd_temperature(mass,vel)*rescale_temp 270 271 ABI_DEALLOCATE(hxredpoint) 272 ABI_DEALLOCATE(xredpoint) 273 ABI_DEALLOCATE(forces_pimd_red) 274 ABI_DEALLOCATE(fsup) 275 276 else 277 278 ! ========================= OTHER TIME STEPS ====================================== 279 280 ! Additional allocations 281 ABI_ALLOCATE(hxredpoint,(3,natom,trotter)) 282 ABI_ALLOCATE(xredpoint,(3,natom,trotter)) 283 ABI_ALLOCATE(forces_pimd_red,(3,natom)) 284 ABI_ALLOCATE(fsup,(3,natom)) 285 286 ! first estimation of ddh, pg and its trace: 287 call matr3inv(rprimd,invrprimd) 288 pg=wg*matmul(ddh,invrprimd) 289 tracepg=pg(1,1)+pg(2,2)+pg(3,3) 290 291 ! Momenta hxredpoint = H ds/dt: estimation 292 if (itimimage==2) then 293 hxredpoint=vel 294 else 295 do iimage=1,trotter 296 hxredpoint(:,:,iimage)=matmul(rprimd,vel(:,:,iimage)) 297 ! because vel in entrance is a scaled velocity (ds/dt) 298 end do 299 end if 300 301 ! Compute temperature at t 302 temperature1=pimd_temperature(mass,hxredpoint)*rescale_temp 303 304 ! Generate random numbers 305 call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce) 306 307 ! Generate random numbers for the barostat 308 call pimd_langevin_random_bar(aleabar,pimd_param%irandom,idum) 309 langev_bar=matmul(aleabar,rprimd)*scalebar 310 311 ! Compute PIMD and Langevin contributions to forces 312 call pimd_forces(forces,natom,spring,0,trotter,xcart) 313 call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint) 314 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 315 & mass,natom,trotter,pimd_param%wtatcon,xcart) 316 317 ! Compute difference between instantaneous stress and imposed stress (barostat) 318 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart) 319 320 diffstress=pimd_diff_stress(stress_pimd,strtarget) 321 322 ! Compute "force" on supercell vectors 323 dstrhh=matmul(diffstress,rprimd) 324 pgdh=matmul(pg,ddh) 325 fg(:,:)=volume*dstrhh(:,:)+pgdh(:,:)+temperature1*kb_HaK*rprimd(:,:) 326 327 ! Evolve the supercell (fist estimation) 328 rprimd_next=two*rprimd-rprimd_prev+(fg-wg*frictionbar*ddh+langev_bar)*dtion*dtion/wg 329 330 ! Evolve atomic positions (first estimation) 331 tmp=matmul(invrprimd,ddh) 332 do iimage=1,trotter 333 call xcart2xred(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage)) 334 fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage)) 335 call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red) 336 337 do iatom=1,natom 338 xred_next(:,iatom,iimage)= & 339 & two*xred(:,iatom,iimage) - xred_prev(:,iatom,iimage) & 340 & +(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom) & 341 & -tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) & 342 & *dtion*dtion/inertmass(iatom) 343 end do 344 end do 345 346 ! Self-consistent loop 347 temperature2=pimd_temperature(mass,xredpoint)*rescale_temp 348 temp2_prev=temperature2; tol=one 349 350 do while (tol>tolerance) 351 ! Reestimate dH/dt at t 352 ddh(:,:)=(rprimd_next(:,:)-rprimd_prev(:,:))/(two*dtion) 353 354 ! Reestimate the scaled velocities at t 355 do iimage=1,trotter 356 xredpoint(:,:,iimage)=(xred_next(:,:,iimage)-xred_prev(:,:,iimage))/(two*dtion) 357 call xred2xcart(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage)) 358 end do 359 ! Reestimate the forces 360 call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,hxredpoint) 361 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 362 & mass,natom,trotter,pimd_param%wtatcon,xcart) 363 ! Compute variation of temperature (to check convergence of SC loop) 364 temperature2=pimd_temperature(mass,xredpoint)*rescale_temp 365 tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev) 366 temp2_prev=temperature2 367 ! Recompute the temperature 368 temperature2=pimd_temperature(mass,hxredpoint)*rescale_temp 369 ! Recompute pg 370 pg=wg*matmul(ddh,invrprimd) 371 tracepg=pg(1,1)+pg(2,2)+pg(3,3) 372 ! Recompute difference between instantaneous stress and imposed stress (barostat) 373 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,hxredpoint,volume,xcart) 374 375 diffstress=pimd_diff_stress(stress_pimd,strtarget) 376 377 ! Recompute "force" on supercell vectors 378 dstrhh=matmul(diffstress,rprimd) 379 pgdh=matmul(pg,ddh) 380 fg(:,:)=volume*diffstress(:,:)+pgdh(:,:)+temperature2*kb_HaK*rprimd(:,:) 381 ! Evolve the supercell (better estimation) 382 rprimd_next=two*rprimd-rprimd_prev+(fg-wg*frictionbar*ddh+langev_bar)*dtion*dtion/wg 383 384 ! Evolve atomic positions (better estimation): 385 tmp=matmul(invrprimd,ddh) 386 do iimage=1,trotter 387 call xcart2xred(natom,rprimd,hxredpoint(:,:,iimage),xredpoint(:,:,iimage)) 388 fsup(:,:)=matmul(tmp,xredpoint(:,:,iimage)) 389 call xcart2xred(natom,rprimd,forces_pimd(:,:,iimage),forces_pimd_red) 390 do iatom=1,natom 391 xred_next(:,iatom,iimage)= & 392 & two*xred(:,iatom,iimage) - xred_prev(:,iatom,iimage) & 393 & +(forces_pimd_red(:,iatom)-two*inertmass(iatom)*fsup(:,iatom) & 394 & -tracepg*inertmass(iatom)*xredpoint(:,iatom,iimage)/(wg*dble(ndof))) & 395 & *dtion*dtion/inertmass(iatom) 396 end do 397 end do 398 end do ! End self-consistent loop 399 400 ! Computation of true temperature from true velocities at t 401 do iimage=1,trotter 402 vel(:,:,iimage)=hxredpoint(:,:,iimage)+matmul(ddh,xred(:,:,iimage)) 403 end do 404 temperature2=pimd_temperature(mass,vel)*rescale_temp 405 406 ! Computation of the real stress tensor at t 407 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart) 408 409 ! translate from pressure to stress 410 stress_pimd=-stress_pimd 411 412 ! Deallocations (Verlet algo) 413 ABI_DEALLOCATE(xredpoint) 414 ABI_DEALLOCATE(hxredpoint) 415 ABI_DEALLOCATE(forces_pimd_red) 416 ABI_DEALLOCATE(fsup) 417 418 end if ! itimimage==1 419 420 !############# Final operations ############################ 421 422 !Compute contributions to energy 423 call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring,trotter,xcart) 424 425 !Print messages 426 prtstress=1 427 call pimd_print(constraint,constraint_output,& 428 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,& 429 & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,& 430 & stress_pimd,temperature1,temperature2,& 431 & pimd_param%traj_unit,trotter,vel,ddh,xcart,xred) 432 433 if (itimimage>1) then 434 ! Estimation of ds/dt at t+dt 435 vel = (three*xred_next - four*xred + xred_prev)/(two * dtion) 436 ddh = (three*rprimd_next - four*rprimd + rprimd_prev)/(two * dtion) 437 end if 438 439 !Come back to reduced coordinates 440 if(itimimage<=1) then 441 do iimage=1,trotter 442 call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage)) 443 end do 444 end if 445 446 !Return cell velocities (does not depend on Trotter) 447 do iimage=1,trotter 448 vel_cell(:,:,iimage)=ddh(:,:) 449 end do 450 451 !Free memory 452 ABI_DEALLOCATE(xcart) 453 ABI_DEALLOCATE(xcart_prev) 454 ABI_DEALLOCATE(xcart_next) 455 ABI_DEALLOCATE(forces_orig) 456 ABI_DEALLOCATE(forces_pimd) 457 ABI_DEALLOCATE(inertmass) 458 ABI_DEALLOCATE(quantummass) 459 ABI_DEALLOCATE(mass) 460 ABI_DEALLOCATE(spring) 461 ABI_DEALLOCATE(alea) 462 ABI_DEALLOCATE(langev) 463 464 end subroutine pimd_langevin_npt