TABLE OF CONTENTS
ABINIT/m_pimd_nosehoover [ Modules ]
NAME
m_pimd_nosehoover
FUNCTION
COPYRIGHT
Copyright (C) 2011-2024 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_nosehoover 23 24 use defs_basis 25 use m_pimd 26 use m_abicore 27 28 use m_geometry, only : xcart2xred, xred2xcart 29 30 implicit none 31 32 private
ABINIT/pimd_nosehoover_npt [ Functions ]
NAME
pimd_nosehoover_npt
FUNCTION
Predicts new positions in Path Integral Molecular Dynamics using Nose-Hoover 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 No change of acell and rprim at present.
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 + Langevin 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
SOURCE
86 subroutine pimd_nosehoover_npt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,& 87 & rprimd,rprimd_next,rprimd_prev,stressin,trotter,vel,vel_cell,& 88 & volume,xred,xred_next,xred_prev) 89 90 implicit none 91 92 !Arguments ------------------------------------ 93 !scalars 94 integer,intent(in) :: itimimage,natom,prtvolimg,trotter 95 real(dp),intent(in) :: volume 96 type(pimd_type),intent(in) :: pimd_param 97 !arrays 98 real(dp),intent(in) :: etotal(trotter),rprimd(3,3),rprimd_prev(3,3),stressin(3,3,trotter) 99 real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter) 100 real(dp),intent(out) :: rprimd_next(3,3),xred_next(3,natom,trotter) 101 real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter),vel_cell(3,3) 102 103 !Local variables------------------------------- 104 !Options 105 real(dp),parameter :: tolerance=tol7 ! SCF tolerance 106 !scalars 107 integer :: idum=-5 108 integer :: constraint,iimage,irestart,ndof,nnos,pitransform,prtstress 109 real(dp) :: dtion,eharm,eharm2,epot,initemp,kt,temperature1,temperature2,thermtemp 110 !arrays 111 real(dp) :: constraint_output(2),ddh(3,3),stress_pimd(3,3,3) 112 real(dp),allocatable :: dzeta(:,:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:) 113 real(dp),allocatable :: inertmass(:),masseff(:,:),qmass(:),quantummass(:),springeff(:,:) 114 real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:),zeta(:) 115 116 ! ************************************************************************* 117 118 !############# Initializations ########################### 119 120 !Allocation of local arrays 121 ABI_MALLOC(xcart,(3,natom,trotter)) 122 ABI_MALLOC(xcart_prev,(3,natom,trotter)) 123 ABI_MALLOC(xcart_next,(3,natom,trotter)) 124 ABI_MALLOC(forces_orig,(3,natom,trotter)) 125 ABI_MALLOC(forces_pimd,(3,natom,trotter)) 126 ABI_MALLOC(inertmass,(natom)) 127 ABI_MALLOC(quantummass,(natom)) 128 129 !Fill in the local variables 130 ndof=3*natom*trotter 131 quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass 132 inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass 133 initemp=pimd_param%mdtemp(1);thermtemp=pimd_param%mdtemp(2) 134 dtion=pimd_param%dtion;pitransform=pimd_param%pitransform 135 kt=thermtemp*kb_HaK 136 forces_orig=forces 137 138 !Allocation/initialization of local variables used for Nose-Hoover chains 139 !Associated variables: 140 !nnos = number of thermostats 141 !dzeta(3,natom,trotter,nnos) = variables of thermostats, in (atomic time unit)^(-1) 142 !qmass(nnos) = masses of thermostats 143 !specific to PIMD: pitransform = coordinate transformation (0:no; 1:normal mode; 2:staging) 144 nnos=pimd_param%nnos 145 ABI_MALLOC(qmass,(nnos)) 146 ABI_MALLOC(zeta,(nnos)) 147 ABI_MALLOC(dzeta,(3,natom,trotter,nnos)) 148 qmass(1:nnos)=pimd_param%qmass(1:nnos) 149 zeta=zero;dzeta=zero 150 151 !Compute cartesian coordinates 152 do iimage=1,trotter 153 call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage)) 154 call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage)) 155 end do 156 157 !Determine if it is a restart or not 158 !If this is a calculation from scratch,generate random distribution of velocities 159 irestart=1;if (itimimage==1) irestart=pimd_is_restart(masseff,vel,vel_cell) 160 161 !Initialize derivatives 162 if (mod(irestart,10)==0) then 163 call pimd_initvel(idum,masseff,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon) 164 end if 165 !vel_cell does not depend on Trotter... 166 ddh=vel_cell(:,:);if (irestart<10) ddh=zero 167 168 !Compute temperature at t 169 temperature1=pimd_temperature(masseff,vel) 170 171 !################## Images evolution ##################### 172 173 !This is temporary 174 xcart_next=zero 175 rprimd_next=rprimd_prev 176 temperature2=pimd_temperature(masseff,vel) 177 178 !Compute contributions to energy 179 call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,springeff,trotter,xcart) 180 181 !Compute stress tensor at t from virial theorem 182 call pimd_stresses(masseff,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart) 183 184 185 !############# Final operations ############################ 186 187 !Print messages 188 prtstress=1 189 call pimd_print(constraint,constraint_output,& 190 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,& 191 & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,& 192 & stress_pimd,temperature1,temperature2,& 193 & pimd_param%traj_unit,trotter,vel,ddh,xcart,xred) 194 195 !Come back to reduced coordinates 196 do iimage=1,trotter 197 call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage)) 198 end do 199 200 !Return cell velocities (does not depend on Trotter) 201 vel_cell(:,:)=ddh(:,:) 202 203 !Free memory 204 ABI_FREE(xcart) 205 ABI_FREE(xcart_prev) 206 ABI_FREE(xcart_next) 207 ABI_FREE(forces_orig) 208 ABI_FREE(forces_pimd) 209 ABI_FREE(inertmass) 210 ABI_FREE(quantummass) 211 ABI_FREE(masseff) 212 ABI_FREE(springeff) 213 ABI_FREE(qmass) 214 ABI_FREE(dzeta) 215 ABI_FREE(zeta) 216 217 end subroutine pimd_nosehoover_npt
ABINIT/pimd_nosehoover_nvt [ Functions ]
NAME
pimd_nosehoover_nvt
FUNCTION
Predicts new positions in Path Integral Molecular Dynamics using Nose-Hoover 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 No change of acell and rprim at present.
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=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
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
Thermization by Nose-Hoover chains according to Martyna, Klein, Tuckerman, J. Chem. Phys. 97, 2635 (1992) [[cite:Martyna1992]] Tuckerman, Marx, Klein, Parrinello, J. Chem. Phys. 104, 5579 (1996) [[cite:Tuckerman1996]]
SOURCE
264 subroutine pimd_nosehoover_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,& 265 & rprimd,stressin,trotter,vel,volume,xred,xred_next,xred_prev) 266 267 implicit none 268 269 !Arguments ------------------------------------ 270 !scalars 271 integer,intent(in) :: itimimage,natom,prtvolimg,trotter 272 real(dp),intent(in) :: volume 273 type(pimd_type),intent(inout) :: pimd_param 274 !arrays 275 real(dp),intent(in) :: etotal(trotter),rprimd(3,3),stressin(3,3,trotter) 276 real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter) 277 real(dp),intent(out) :: xred_next(3,natom,trotter) 278 real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter) 279 280 !Local variables------------------------------- 281 !Options 282 real(dp),parameter :: tolerance=tol9 ! SCF tolerance 283 !scalars 284 integer :: idum=-5 285 integer :: iimage,irestart,ndof,nnos,pitransform,prtstress 286 real(dp) :: dtion,eharm,eharm2,epot,initemp,kt 287 real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol 288 character(len=500) :: msg 289 !arrays 290 real(dp) :: constraint_output(2),spring_prim(natom),stress_pimd(3,3,3),vel_cell(3,3) 291 real(dp),allocatable :: forces_orig(:,:,:),forces_pimd(:,:,:) 292 real(dp),allocatable :: inertmass(:),mass(:,:),qmass(:),quantummass(:),spring(:,:) 293 real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:) 294 real(dp),allocatable :: dzeta(:,:,:,:),zeta_prev(:,:,:,:),zeta(:,:,:,:) 295 real(dp),allocatable :: zeta_next(:,:,:,:) 296 297 ! ************************************************************************* 298 299 !############# Initializations ########################### 300 301 !Allocation of local arrays 302 ABI_MALLOC(xcart,(3,natom,trotter)) 303 ABI_MALLOC(xcart_prev,(3,natom,trotter)) 304 ABI_MALLOC(xcart_next,(3,natom,trotter)) 305 ABI_MALLOC(forces_orig,(3,natom,trotter)) 306 ABI_MALLOC(forces_pimd,(3,natom,trotter)) 307 ABI_MALLOC(inertmass,(natom)) 308 ABI_MALLOC(quantummass,(natom)) 309 310 !Fill in the local variables 311 ndof=3*natom*trotter 312 pitransform=pimd_param%pitransform 313 quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass 314 inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass 315 if(pitransform==1) inertmass=quantummass !compulsory for good definition of normal mode masses 316 if(pitransform==2) inertmass=quantummass !compulsory for good definition of staging masses 317 initemp=pimd_param%mdtemp(1);thermtemp=pimd_param%mdtemp(2) 318 dtion=pimd_param%dtion 319 kt=thermtemp*kb_HaK 320 forces_orig=forces 321 322 !Allocation/initialization of local variables used for Nose-Hoover chains 323 !Associated variables: 324 !nnos = number of thermostats 325 !zeta,zeta_next,zeta_prev(3,natom,trotter,nnos) = variables of thermostats, dzeta, its time derivative 326 !qmass(nnos) = masses of thermostats 327 !specific to PIMD: pitransform = coordinate transformation (0:no; 1:normal mode; 2:staging) 328 nnos=pimd_param%nnos 329 ABI_MALLOC(qmass,(nnos)) 330 ABI_MALLOC(zeta_prev,(3,natom,trotter,nnos)) 331 ABI_MALLOC(zeta,(3,natom,trotter,nnos)) 332 ABI_MALLOC(zeta_next,(3,natom,trotter,nnos)) 333 ABI_MALLOC(dzeta,(3,natom,trotter,nnos)) 334 !initialization 335 qmass(1:nnos)=pimd_param%qmass(1:nnos) 336 zeta_prev(:,:,:,:)=pimd_param%zeta_prev(:,:,:,:) 337 zeta(:,:,:,:) =pimd_param%zeta(:,:,:,:) 338 dzeta(:,:,:,:) =pimd_param%dzeta(:,:,:,:) !unuseful to initialize zeta_next 339 340 !Masses and spring constants (according to pitransform) 341 select case(pitransform) 342 case(0) 343 ABI_MALLOC(mass,(natom,1)) 344 ABI_MALLOC(spring,(natom,1)) 345 case(1,2) 346 ABI_MALLOC(mass,(natom,trotter)) 347 ABI_MALLOC(spring,(natom,trotter)) 348 end select 349 spring_prim(:)=quantummass(:)*dble(trotter)*kt*kt 350 351 call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,pitransform,trotter) 352 353 !Recommended value of Nose mass 354 write(msg,'(2a,f9.2,3a)') ch10,& 355 & ' Recommended value of Nose mass is',one/(dble(trotter)*kt),' (atomic units)',ch10,& 356 & '(see Tuckerman et al, J. Chem. Phys. 104, 5579 (1996))' ! [[cite:Tuckerman1996]] 357 call wrtout(std_out,msg,'COLL') 358 359 !Compute cartesian coordinates 360 do iimage=1,trotter 361 call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage)) 362 call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage)) 363 end do 364 365 !Determine if it is a restart or not 366 !If this is a calculation from scratch,generate random distribution of velocities 367 irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel) 368 if (irestart==0) then 369 call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon) 370 end if 371 372 !Compute temperature at t 373 temperature1=pimd_temperature(mass,vel) 374 375 !################## Images evolution ##################### 376 377 !Transform the coordinates and forces (according to pitransform) 378 call pimd_coord_transform(xcart,1,natom,pitransform,trotter) 379 call pimd_force_transform(forces,1,natom,pitransform,trotter) !compute staging forces 380 call pimd_forces(forces,natom,spring,pitransform,trotter,xcart) 381 call pimd_nosehoover_forces(dzeta,forces,forces_pimd,mass,natom,nnos,trotter,vel) 382 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 383 & mass,natom,trotter,pimd_param%wtatcon,xcart) 384 385 !Compute atomic positions at t+dt 386 if (itimimage<=1) then 387 388 ! === 1st time step: single Taylor algorithm 389 ! Predict positions 390 call pimd_predict_taylor(dtion,forces_pimd,mass,natom,trotter,& 391 & vel,xcart,xcart_next) 392 393 ! Estimate the velocities at t+dt/2 394 vel=(xcart_next-xcart)/dtion 395 396 ! Compute new temperature 397 temperature2=pimd_temperature(mass,vel) 398 399 ! Propagate the thermostat variables 400 call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,& 401 & thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform) 402 403 dzeta=(zeta_next-zeta)/dtion 404 405 else 406 407 ! === Other time steps: Verlet algorithm + SC cycle 408 ! Predict positions 409 call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter) 410 call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,& 411 & xcart,xcart_next,xcart_prev) 412 ! Propagate the thermostat variables 413 call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,& 414 & thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform) 415 ! Self-consistent loop 416 temperature2=pimd_temperature(mass,vel) 417 temp2_prev=temperature2; tol=one 418 do while (tol>tolerance) 419 ! Recompute a (better) estimation of the velocity at time step t 420 vel = (xcart_next - xcart_prev) / (two*dtion) 421 dzeta=(zeta_next - zeta_prev) / (two*dtion) 422 temperature2=pimd_temperature(mass,vel) 423 ! Reestimate the force 424 call pimd_nosehoover_forces(dzeta,forces,forces_pimd,mass,natom,nnos,trotter,vel) 425 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 426 & mass,natom,trotter,pimd_param%wtatcon,xcart) 427 ! Compute new positions 428 call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,& 429 & xcart,xcart_next,xcart_prev) 430 ! Propagate the thermostat variables 431 call pimd_nosehoover_propagate(dtion,dzeta,mass,natom,nnos,qmass,& 432 & thermtemp,trotter,vel,zeta,zeta_next,zeta_prev,itimimage,pitransform) 433 ! Compute variation of temperature (to check convergence of SC loop) 434 tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev) 435 temp2_prev=temperature2 436 end do ! End self-consistent loop 437 438 end if ! itimimage==1 439 440 !Come back to primitive coordinates and velocities 441 call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter) 442 call pimd_coord_transform(xcart ,-1,natom,pitransform,trotter) 443 call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter) 444 445 !Compute contributions to energy 446 call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring_prim,trotter,xcart) 447 448 !Compute stress tensor at t from virial theorem 449 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart) 450 stress_pimd=-stress_pimd ! Translate pressure to stress 451 452 !############# Final operations ############################ 453 454 !Print messages 455 vel_cell=zero;prtstress=1;if (prtvolimg>=2) prtstress=0 456 call pimd_print(pimd_param%constraint,constraint_output,& 457 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,& 458 & itimimage,kt,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,& 459 & stress_pimd,temperature1,temperature2,& 460 & pimd_param%traj_unit,trotter,vel,vel_cell,xcart,xred) 461 462 !If possible, estimate the (transformed) velocities at t+dt 463 if (itimimage>1) then 464 call pimd_coord_transform(xcart_next,1,natom,pitransform,trotter) 465 call pimd_coord_transform(xcart,1,natom,pitransform,trotter) 466 call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter) 467 vel = (three*xcart_next - four*xcart + xcart_prev)/(two * dtion) 468 call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter) 469 call pimd_coord_transform(xcart,-1,natom,pitransform,trotter) 470 call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter) 471 end if 472 473 !Come back to reduced coordinates 474 do iimage=1,trotter 475 call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage)) 476 end do 477 478 !update thermostat variables 479 dzeta = (three*zeta_next - four*zeta + zeta_prev)/(two * dtion) 480 zeta_prev=zeta 481 zeta=zeta_next 482 pimd_param%zeta_prev(:,:,:,:)=zeta_prev(:,:,:,:) 483 pimd_param%zeta(:,:,:,:) =zeta(:,:,:,:) 484 pimd_param%dzeta(:,:,:,:) =dzeta(:,:,:,:) 485 486 !Free memory 487 ABI_FREE(xcart) 488 ABI_FREE(xcart_prev) 489 ABI_FREE(xcart_next) 490 ABI_FREE(forces_orig) 491 ABI_FREE(forces_pimd) 492 ABI_FREE(inertmass) 493 ABI_FREE(quantummass) 494 ABI_FREE(mass) 495 ABI_FREE(spring) 496 ABI_FREE(qmass) 497 ABI_FREE(zeta_prev) 498 ABI_FREE(zeta) 499 ABI_FREE(zeta_next) 500 ABI_FREE(dzeta) 501 502 end subroutine pimd_nosehoover_nvt