1from ase.units import Bohr, AUT, _me, _amu 2from gpaw.tddft.units import attosec_to_autime 3from gpaw.forces import calculate_forces 4 5############################################################################### 6# EHRENFEST DYNAMICS WITHIN THE PAW METHOD 7# WORKS WITH PAW AS LONG THERE ISN'T TOO 8# TOO MUCH OVERLAP BETWEEN THE SPHERES 9# SUPPORTS ALSO HGH PSEUDOPOTENTIALS 10############################################################################### 11 12# m a(t+dt) = F[psi(t),x(t)] 13# x(t+dt/2) = x(t) + v(t) dt/2 + .5 a(t) (dt/2)^2 14# vh(t+dt/2) = v(t) + .5 a(t) dt/2 15# m a(t+dt/2) = F[psi(t),x(t+dt/2)] 16# v(t+dt/2) = vh(t+dt/2) + .5 a(t+dt/2) dt/2 17# 18# psi(t+dt) = U(t,t+dt) psi(t) 19# 20# m a(t+dt/2) = F[psi(t+dt),x(t+dt/2)] 21# x(t+dt) = x(t+dt/2) + v(t+dt/2) dt/2 + .5 a(t+dt/2) (dt/2)^2 22# vh(t+dt) = v(t+dt/2) + .5 a(t+dt/2) dt/2 23# m a(t+dt) = F[psi(t+dt),x(t+dt)] 24# v(t+dt) = vh(t+dt/2) + .5 a(t+dt/2) dt/2 25 26# TODO: move force corrections from forces.py to this module, as well as 27# the cg method for calculating the inverse of S from overlap.py 28 29 30class EhrenfestVelocityVerlet: 31 32 def __init__(self, calc, mass_scale=1.0, setups='paw'): 33 """Initializes the Ehrenfest MD calculator. 34 35 Parameters 36 ---------- 37 38 calc: TDDFT Object 39 40 mass_scale: 1.0 41 Scaling coefficient for atomic masses 42 43 setups: {'paw', 'hgh'} 44 Type of setups to use for the calculation 45 46 Note 47 ------ 48 49 Use propagator = 'EFSICN' for when creating the TDDFT object from a 50 PAW ground state calculator and propagator = 'EFSICN_HGH' for HGH 51 pseudopotentials 52 53 """ 54 self.calc = calc 55 self.setups = setups 56 self.x = self.calc.atoms.positions.copy() / Bohr 57 self.xn = self.x.copy() 58 self.v = self.x.copy() 59 amu_to_aumass = _amu / _me 60 if self.calc.atoms.get_velocities() is not None: 61 self.v = self.calc.atoms.get_velocities() / (Bohr / AUT) 62 else: 63 self.v[:] = 0.0 64 self.calc.atoms.set_velocities(self.v) 65 66 self.vt = self.v.copy() 67 self.vh = self.v.copy() 68 self.time = 0.0 69 70 self.M = calc.atoms.get_masses() * amu_to_aumass * mass_scale 71 72 self.a = self.v.copy() 73 self.ah = self.a.copy() 74 self.an = self.a.copy() 75 self.F = self.a.copy() 76 77 self.calc.get_td_energy() 78 self.F = self.get_forces() 79 80 for i in range(len(self.F)): 81 self.a[i] = self.F[i] / self.M[i] 82 83 def get_forces(self): 84 return calculate_forces(self.calc.wfs, 85 self.calc.td_density.get_density(), 86 self.calc.td_hamiltonian.hamiltonian, 87 self.calc.log) 88 89 def propagate(self, dt): 90 """Performs one Ehrenfest MD propagation step 91 92 Parameters 93 ---------- 94 95 dt: scalar 96 Time step (in attoseconds) used for the Ehrenfest MD step 97 98 """ 99 self.x = self.calc.atoms.positions.copy() / Bohr 100 self.v = self.calc.atoms.get_velocities() / (Bohr / AUT) 101 102 dt = dt * attosec_to_autime 103 104 # m a(t+dt) = F[psi(t),x(t)] 105 self.calc.atoms.positions = self.x * Bohr 106 self.calc.set_positions(self.calc.atoms) 107 self.calc.get_td_energy() 108 self.F = self.get_forces() 109 110 for i in range(len(self.F)): 111 self.a[i] = self.F[i] / self.M[i] 112 113 # x(t+dt/2) = x(t) + v(t) dt/2 + .5 a(t) (dt/2)^2 114 # vh(t+dt/2) = v(t) + .5 a(t) dt/2 115 self.xh = self.x + self.v * dt / 2 + .5 * self.a * dt / 2 * dt / 2 116 self.vhh = self.v + .5 * self.a * dt / 2 117 118 # m a(t+dt/2) = F[psi(t),x(t+dt/2)a] 119 self.calc.atoms.positions = self.xh * Bohr 120 self.calc.set_positions(self.calc.atoms) 121 self.calc.get_td_energy() 122 self.F = self.get_forces() 123 124 for i in range(len(self.F)): 125 self.ah[i] = self.F[i] / self.M[i] 126 127 # v(t+dt/2) = vh(t+dt/2) + .5 a(t+dt/2) dt/2 128 self.vh = self.vhh + 0.5 * self.ah * dt / 2 129 130 # Propagate wf 131 # psi(t+dt) = U(t,t+dt) psi(t) 132 self.propagate_single(dt) 133 134 # m a(t+dt/2) = F[psi(t+dt),x(t+dt/2)] 135 self.calc.atoms.positions = self.xh * Bohr 136 self.calc.set_positions(self.calc.atoms) 137 self.calc.get_td_energy() 138 self.F = self.get_forces() 139 140 for i in range(len(self.F)): 141 self.ah[i] = self.F[i] / self.M[i] 142 143 # x(t+dt) = x(t+dt/2) + v(t+dt/2) dt/2 + .5 a(t+dt/2) (dt/2)^2 144 # vh(t+dt) = v(t+dt/2) + .5 a(t+dt/2) dt/2 145 self.xn = self.xh + self.vh * dt / 2 + 0.5 * self.ah * dt / 2 * dt / 2 146 self.vhh = self.vh + .5 * self.ah * dt / 2 147 148 # m a(t+dt) = F[psi(t+dt),x(t+dt)] 149 self.calc.atoms.positions = self.xn * Bohr 150 self.calc.set_positions(self.calc.atoms) 151 self.calc.get_td_energy() 152 self.calc.update_eigenvalues() 153 self.F = self.get_forces() 154 155 for i in range(len(self.F)): 156 self.an[i] = self.F[i] / self.M[i] 157 158 # v(t+dt) = vh(t+dt/2) + .5 a(t+dt/2) dt/2 159 self.vn = self.vhh + .5 * self.an * dt / 2 160 161 # update 162 self.x[:] = self.xn 163 self.v[:] = self.vn 164 self.a[:] = self.an 165 166 # update atoms 167 self.calc.atoms.set_positions(self.x * Bohr) 168 self.calc.atoms.set_velocities(self.v * Bohr / AUT) 169 170 def propagate_single(self, dt): 171 if self.setups == 'paw': 172 self.calc.propagator.propagate(self.time, dt, self.vh) 173 else: 174 self.calc.propagator.propagate(self.time, dt) 175 176 def get_energy(self): 177 """Updates kinetic, electronic and total energies""" 178 179 self.Ekin = 0.5 * (self.M * (self.v**2).sum(axis=1)).sum() 180 self.e_coulomb = self.calc.get_td_energy() 181 self.Etot = self.Ekin + self.e_coulomb 182 return self.Etot 183 184 def get_velocities_in_au(self): 185 return self.v 186 187 def set_velocities_in_au(self, v): 188 self.v[:] = v 189 self.calc.atoms.set_velocities(v * Bohr / AUT) 190