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