TABLE OF CONTENTS
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
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) 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
See Quigley,Probert, JCP 120, 11432 (2004), part III
PARENTS
predict_pimd
CHILDREN
pimd_apply_constraint,pimd_coord_transform,pimd_energies pimd_force_transform,pimd_forces,pimd_initvel,pimd_langevin_forces pimd_langevin_random,pimd_langevin_random_init,pimd_langevin_random_qtb pimd_mass_spring,pimd_predict_taylor,pimd_predict_verlet,pimd_print pimd_stresses,xcart2xred,xred2xcart
SOURCE
61 #if defined HAVE_CONFIG_H 62 #include "config.h" 63 #endif 64 65 #include "abi_common.h" 66 67 subroutine pimd_langevin_nvt(etotal,forces,itimimage,natom,pimd_param,prtvolimg,& 68 & rprimd,stressin,trotter,vel,volume,xred,xred_next,xred_prev) 69 70 use defs_basis 71 use m_errors 72 use m_pimd 73 use m_profiling_abi 74 75 use m_geometry, only : xcart2xred, xred2xcart 76 77 !This section has been created automatically by the script Abilint (TD). 78 !Do not modify the following lines by hand. 79 #undef ABI_FUNC 80 #define ABI_FUNC 'pimd_langevin_nvt' 81 !End of the abilint section 82 83 implicit none 84 85 !Arguments ------------------------------------ 86 !scalars 87 integer,intent(in) :: itimimage,natom,prtvolimg,trotter 88 real(dp),intent(in) :: volume 89 type(pimd_type),intent(in) :: pimd_param 90 !arrays 91 real(dp),intent(in) :: etotal(trotter),rprimd(3,3),stressin(3,3,trotter) 92 real(dp),intent(in),target :: xred(3,natom,trotter),xred_prev(3,natom,trotter) 93 real(dp),intent(out) :: xred_next(3,natom,trotter) 94 real(dp),intent(inout) :: forces(3,natom,trotter),vel(3,natom,trotter) 95 96 !Local variables------------------------------- 97 !Options 98 ! Option for the Langevin algorithm correction 99 integer,parameter :: ilangevin=0 100 ! Tolerance for the SC cycle 101 real(dp),parameter :: tolerance=tol9 102 103 !scalars 104 integer :: idum=-5 105 integer :: iimage,irestart,ndof,pitransform,prtstress,use_qtb,zeroforce 106 real(dp) :: dtion,eharm,eharm2,epot,friction,initemp,kt,kt_,rescale_temp 107 real(dp) :: temperature1,temperature2,temp2_prev,thermtemp,tol 108 !arrays 109 real(dp) :: constraint_output(2),spring_prim(natom),stress_pimd(3,3,3),vel_cell(3,3) 110 real(dp),allocatable :: alea(:,:,:),forces_orig(:,:,:),forces_pimd(:,:,:) 111 real(dp),allocatable :: inertmass(:),langev(:,:),mass(:,:),quantummass(:),spring(:,:) 112 real(dp),allocatable :: xcart(:,:,:),xcart_next(:,:,:),xcart_prev(:,:,:) 113 114 ! ************************************************************************* 115 116 if (pimd_param%use_qtb==1.and.pimd_param%qtb_file_unit<=0) then 117 MSG_BUG('piqtb_force not open!') 118 end if 119 120 !############# Initializations ########################### 121 122 pitransform=pimd_param%pitransform 123 124 !The following option forces the total of forces to be zero 125 !It prevents the translation of the center of mass 126 !If it is zero, no constraint on mass center is applied 127 zeroforce=1 128 if(pitransform==1) zeroforce=0 129 if(pitransform==2) zeroforce=0 130 if(pimd_param%constraint==1) zeroforce=0 131 132 !Allocation of local arrays 133 ABI_ALLOCATE(xcart,(3,natom,trotter)) 134 ABI_ALLOCATE(xcart_prev,(3,natom,trotter)) 135 ABI_ALLOCATE(xcart_next,(3,natom,trotter)) 136 ABI_ALLOCATE(forces_orig,(3,natom,trotter)) 137 ABI_ALLOCATE(forces_pimd,(3,natom,trotter)) 138 ABI_ALLOCATE(inertmass,(natom)) 139 ABI_ALLOCATE(quantummass,(natom)) 140 141 !Fill in the local variables 142 use_qtb=pimd_param%use_qtb 143 ndof=3*natom*trotter 144 rescale_temp=one;if(zeroforce==1)rescale_temp=dble(ndof)/dble(ndof-3) 145 quantummass(1:natom)=pimd_param%amu (pimd_param%typat(1:natom))*amu_emass 146 inertmass (1:natom)=pimd_param%pimass(pimd_param%typat(1:natom))*amu_emass 147 if(pitransform==1) inertmass=quantummass !compulsory for good definition of normal mode masses 148 if(pitransform==2) inertmass=quantummass !compulsory for good definition of staging masses 149 initemp=pimd_param%mdtemp(1)/rescale_temp;thermtemp=pimd_param%mdtemp(2) 150 friction=pimd_param%vis;dtion=pimd_param%dtion 151 kt=thermtemp*kb_HaK 152 forces_orig=forces 153 154 !Masses and spring constants 155 select case(pitransform) 156 case(0) 157 ABI_ALLOCATE(mass,(natom,1)) ! This second dimension is needed 158 ABI_ALLOCATE(spring,(natom,1)) 159 ABI_ALLOCATE(langev,(natom,1)) 160 case(1,2) 161 ABI_ALLOCATE(mass,(natom,trotter)) 162 ABI_ALLOCATE(spring,(natom,trotter)) 163 ABI_ALLOCATE(langev,(natom,trotter)) 164 end select 165 spring_prim(:)=quantummass(:)*dble(trotter)*kt*kt 166 call pimd_mass_spring(inertmass,kt,mass,natom,quantummass,spring,pitransform,trotter) 167 168 !Initialize random forces 169 ABI_ALLOCATE(alea,(3,natom,trotter)) 170 if (use_qtb==0) then 171 langev(:,:)=sqrt(two*friction*mass(:,:)*kt/dtion) 172 else 173 langev(:,:)=sqrt(two*friction*mass(:,:)) 174 end if 175 176 !Random number generator initialization 177 if(itimimage<=1) then 178 call pimd_langevin_random_init(pimd_param%irandom,idum) 179 end if 180 181 !Compute cartesian coordinates 182 do iimage=1,trotter 183 call xred2xcart(natom,rprimd,xcart (:,:,iimage),xred(:,:,iimage)) 184 call xred2xcart(natom,rprimd,xcart_prev(:,:,iimage),xred_prev(:,:,iimage)) 185 end do 186 187 !Determine if it is a restart or not 188 !If this is a calculation from scratch,generate random distribution of velocities 189 irestart=1;if (itimimage==1) irestart=pimd_is_restart(mass,vel) 190 if (irestart==0) then 191 call pimd_initvel(idum,mass,natom,initemp,trotter,vel,pimd_param%constraint,pimd_param%wtatcon) 192 end if 193 194 !Compute temperature at t 195 temperature1=pimd_temperature(mass,vel)*rescale_temp 196 197 !################## Images evolution ##################### 198 199 !Generate random numbers 200 if (use_qtb==0) then 201 call pimd_langevin_random(alea,pimd_param%irandom,idum,langev,mass,natom,trotter,zeroforce) 202 else 203 call pimd_langevin_random_qtb(alea,langev,mass,natom,pimd_param%qtb_file_unit,trotter,zeroforce) 204 end if 205 206 !Compute PIMD and Langevin contributions to forces 207 call pimd_coord_transform(xcart,1,natom,pitransform,trotter) 208 call pimd_force_transform(forces,1,natom,pitransform,trotter) !compute staging forces 209 call pimd_forces(forces,natom,spring,pitransform,trotter,xcart) 210 call pimd_langevin_forces(alea,forces,forces_pimd,friction,langev,mass,natom,trotter,vel) 211 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 212 & mass,natom,trotter,pimd_param%wtatcon,xcart) 213 214 !Compute atomic positions at t+dt 215 if (itimimage<=1) then 216 217 ! === 1st time step: single Taylor algorithm 218 ! Predict positions 219 call pimd_predict_taylor(dtion,forces_pimd,mass,natom,trotter,& 220 & vel,xcart,xcart_next) 221 222 ! Estimate the velocities at t+dt/2 223 vel=(xcart_next-xcart)/dtion 224 225 ! Compute new temperature 226 temperature2=pimd_temperature(mass,vel)*rescale_temp 227 228 else 229 230 ! === Other time steps: Verlet algorithm + SC cycle 231 ! Predict positions 232 call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter) 233 call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,& 234 & xcart,xcart_next,xcart_prev) 235 ! Self-consistent loop 236 temperature2=pimd_temperature(mass,vel)*rescale_temp 237 temp2_prev=temperature2; tol=one 238 do while (tol>tolerance) 239 ! Recompute a (better) estimation of the velocity at time step t 240 vel = (xcart_next - xcart_prev) / (two*dtion) 241 temperature2=pimd_temperature(mass,vel)*rescale_temp 242 ! Reestimate the force 243 call pimd_langevin_forces(alea,forces,forces_pimd,friction,& 244 & langev,mass,natom,trotter,vel) 245 call pimd_apply_constraint(pimd_param%constraint,constraint_output,forces_pimd,& 246 & mass,natom,trotter,pimd_param%wtatcon,xcart) 247 ! Compute new positions 248 call pimd_predict_verlet(dtion,forces_pimd,mass,natom,trotter,& 249 & xcart,xcart_next,xcart_prev) 250 251 ! Compute variation of temperature (to check convergence of SC loop) 252 tol=dabs(temperature2-temp2_prev)/dabs(temp2_prev) 253 temp2_prev=temperature2 254 255 end do ! End self-consistent loop 256 257 end if ! itimimage==1 258 259 call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter) 260 call pimd_coord_transform(xcart,-1,natom,pitransform,trotter) 261 call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter) 262 263 !Compute contributions to energy 264 call pimd_energies(eharm,eharm2,epot,etotal,forces_orig,natom,spring_prim,trotter,xcart) 265 266 !Compute stress tensor at t 267 if (use_qtb==0) then 268 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,thermtemp,thermtemp,trotter,vel,volume,xcart) 269 else 270 call pimd_stresses(mass,natom,quantummass,stress_pimd,stressin,temperature2,thermtemp,trotter,vel,volume,xcart) 271 end if 272 stress_pimd=-stress_pimd ! Translate pressure to stress 273 274 !############# Final operations ############################ 275 276 !Print messages 277 vel_cell=zero;prtstress=1;if (prtvolimg>=2) prtstress=0 278 kt_=kt;if (use_qtb==1) kt_=temperature2*kb_HaK 279 call pimd_print(pimd_param%constraint,constraint_output,& 280 & eharm,eharm2,epot,forces_pimd,inertmass,irestart,& 281 & itimimage,kt_,natom,pimd_param%optcell,prtstress,prtvolimg,rprimd,& 282 & stress_pimd,temperature1,temperature2,& 283 & pimd_param%traj_unit,trotter,vel,vel_cell,xcart,xred) 284 285 !If possible, estimate the (transformed) velocities at t+dt 286 if (itimimage>1) then 287 call pimd_coord_transform(xcart_next,1,natom,pitransform,trotter) 288 call pimd_coord_transform(xcart,1,natom,pitransform,trotter) 289 call pimd_coord_transform(xcart_prev,1,natom,pitransform,trotter) 290 vel = (three*xcart_next - four*xcart + xcart_prev)/(two * dtion) 291 call pimd_coord_transform(xcart_next,-1,natom,pitransform,trotter) 292 call pimd_coord_transform(xcart,-1,natom,pitransform,trotter) 293 call pimd_coord_transform(xcart_prev,-1,natom,pitransform,trotter) 294 end if 295 296 !Come back to reduced coordinates 297 do iimage=1,trotter 298 call xcart2xred(natom,rprimd,xcart_next(:,:,iimage),xred_next(:,:,iimage)) 299 end do 300 301 !Free memory 302 ABI_DEALLOCATE(xcart) 303 ABI_DEALLOCATE(xcart_prev) 304 ABI_DEALLOCATE(xcart_next) 305 ABI_DEALLOCATE(forces_orig) 306 ABI_DEALLOCATE(forces_pimd) 307 ABI_DEALLOCATE(inertmass) 308 ABI_DEALLOCATE(quantummass) 309 ABI_DEALLOCATE(mass) 310 ABI_DEALLOCATE(spring) 311 ABI_DEALLOCATE(alea) 312 ABI_DEALLOCATE(langev) 313 314 end subroutine pimd_langevin_nvt