DFT - Solution of the KS equation. trial n(r) build V ... CPMD Car-Parrinello Molecular Dynamics. Stop ... adiabatically the ions remaining close to the ground state. ... with given tempw, press ... ekin_conv_thr = 1.d-4, ... sd -> steepest descend.
Car-Parrinello Molecular Dynamics Simulations With QE Carlo Cavazzoni CINECA Italy
CINECA is a non profit Consortium, made up of 40 Italian universities, The National Institute of Oceanography and Experimental Geophysics - OGS, the CNR (National Research Council), and the Ministry of University and Research.
CINECA is the largest Italian Supercomputing center
Outline CP theory CP basic simulation parameters CP in QE Setup for NVE Simulation Setup for NVT Simulation Setup for NPT Simulation
MD - Molecular Dynamics Initial configuration
accumulate statistics (RDF, Energies, MSD, ecc..)
t=t+dt
Stop
calculate forces update atomic positions
See also http://en.wikipedia.org/wiki/Molecular_dynamics
DFT - Solution of the KS equation. trial n(r)
build VKS
solve [-1/2∇2+VKS]ϕi = εiϕi
nold = nnew
calculate nnew(r) = ∑fi|ϕi|2
no
|nnew-nold| common values 2-10 in Hartree atomic units 1 a.u.=4.8378 * 10^-17 s : beware, PW code use Rydberg atomic units, twice that much!!!)
µ
effective electronic mass
(emass -> common values 100-1000 a.u. of mass 1 a.u. = 1/1822.9 a.m.u. = 9.10939 * 10^-31 kg )
Key CP parameters (for NVT) T
temperature
(tempw -> 300-1000 in Kelvin )
Q
Nose-Hoover thermostat mass
(for the input use the frequency fnosep -> 10-100 in terahertz )
Key CP parameters (for NPT) W
effective cell mass
(wmass ~3ΣMi/4π2 a.u. of mass )
P
external (hydrostatic) pressure
(press -> 0 – 1000 Kbar )
Parameters setup You want to simulate your system with given tempw, press Remaining on the BO surface! Tune: dt, emass, wmass Some test run may be required!
Parameters relations gap
emass
cut-off
dt
wmass
dt
tempw
dt
dt
CP in QE CPV subdirectory Executable: cp.x NO K-points, only GAMMA, no PAW Example: 18, 19, 20, 21, 23, 26, 27, 30, autopilot-example, Restart_example, QExml_example Tests: cptests
A CP simulation require usually many RUNs and JOBs 1. Minimize the Electronic degrees of freedom ( d.o.f. ) 2. Minimize the ionic d.o.f. 3. Randomize ionic d.o.f. 4. Re-minimize Electronic d.o.f 5. Move Electronic and Ionic d.o.f. using Verlet to integrate the equations of motion 6. Change temperature using thermostat 7. Accumulate statistics for few picoseconds of simulated time ( 20000 – 40000 time steps )
Electronic Minimization I/II CP time step (in a.u.) convergence threshold, effective only for minimization
&CONTROL title = ' Water Molecule ', calculation = 'cp', restart_mode = 'from_scratch', prefix = 'h2o_mol‘ nstep = 50, dt = 5.0d0, etot_conv_thr = 1.d-9, ekin_conv_thr = 1.d-4, /
&SYSTEM
Energy cut-off, in Rydberg
ibrav = 1, celldm(1) = 10.0, nat = 3, ntyp = 2, ecutwfc = 70.0, /
number of CP step
Simulation Cell (in a.u.)
number of atom, and species
Electronic Minimization II/II Fictitious electron mass (a.u.) µ parameter in CP dynamic
PW Energy cut-off (in Rydberg) for Fourier Acceleration
&ELECTRONS emass = 400.d0, emass_cutoff = 2.5d0, orthogonalization = 'ortho', electron_dynamics = 'sd', /
electron dynamics sd -> steepest descend damp -> damped verlet -> Verlet
&IONS
ion_dynamics = 'none', /
Label
Mass (uma) PP
Label
x y z
ATOMIC_SPECIES O 16.0d0 O.BLYP.UPF H 1.00d0 H.fpmd.UPF ATOMIC_POSITIONS (bohr) O 0.0099 0.0099 H 1.8325 -0.2243 H -0.2243 1.8325
orthogonalization algorithm ‘ortho’ or 'Gram-Schmidt'
0.0000 -0.0001 0.0002
ion dynamics none -> kept fixed sd -> steepest descend damp -> damped verlet -> Verlet
CP – Dynamics I/II &CONTROL title = ' Water Molecule ', calculation = 'cp', restart_mode = ‘restart', nstep = 50, dt = 5.0d0, prefix = 'h2o_mol' /
verlet should be used for ions and electrons
&IONS
ion_dynamics = ‘verlet', ion_velocities = ‘zero’ /
&SYSTEM ibrav = 1, celldm(1) = 10.0, nat = 3, ntyp = 2, ecutwfc = 70.0, /
ATOMIC_SPECIES O 16.0d0 O.BLYP.UPF 4 H 1.00d0 H.fpmd.UPF 4 ATOMIC_POSITIONS (bohr) O 0.0099 0.0099 H 1.8325 -0.2243 H -0.2243 1.8325
&ELECTRONS emass = 400.d0, emass_cutoff = 2.5d0, orthogonalization = 'ortho', electron_dynamics = ‘verlet', electron_velocities = ‘zero’ /
in the very first run of dynamics velocities should be set to zero
0.0000 -0.0001 0.0002
CP – Dynamics II/II ( thermostat ) &CONTROL title = ' Water Molecule ', calculation = 'cp', restart_mode = ‘restart', nstep = 50, dt = 5.0d0, prefix = 'h2o_mol' / &SYSTEM ibrav = 1, celldm(1) = 10.0, nat = 3, ntyp = 2, ecutwfc = 70.0, / &ELECTRONS emass = 400.d0, emass_cutoff = 2.5d0, orthogonalization = 'ortho', electron_dynamics = ‘verlet', /
Termostat temperature and frequency tempw ( Kelvin ) fnosep ( THz ) &IONS
ion_dynamics = ‘verlet', ion_temperature = ‘nose’ tempw = 300, fnosep = 70.0 / ATOMIC_SPECIES O 16.0d0 O.BLYP.UPF 4 H 1.00d0 H.fpmd.UPF 4 ATOMIC_POSITIONS (bohr) O 0.0099 0.0099 H 1.8325 -0.2243 H -0.2243 1.8325
0.0000 -0.0001 0.0002
Parrinello-Rahman Dynamics &CONTROL title = ' Water Molecule ', calculation = ‘vc-cp', restart_mode = ‘restart', nstep = 50, dt = 5.0d0, prefix = 'h2o_mol' / &SYSTEM ibrav = 1, celldm(1) = 10.0, nat = 3, ntyp = 2, ecutwfc = 80.0, ecfixed = 70.0, qcutz = 70.0, q2sigma = 5.0,
change the calculation type (variable cell cp) &IONS
ion_dynamics = ‘verlet', Add the CELL namelists, with cell dynamics: pr -> parrinello rahman sd -> steepest descent damp -> damped Costant cutoff issue
/ &ELECTRONS emass = 400.d0, emass_cutoff = 2.5d0, orthogonalization = 'ortho', electron_dynamics = ‘verlet', /
external pressure (Kbar)
/ &CELL
cell_dynamics = ‘pr’, press = 0.0 / ATOMIC_SPECIES O 16.0d0 O.BLYP.UPF 4 H 1.00d0 H.fpmd.UPF 4 ATOMIC_POSITIONS (bohr) O 0.0099 0.0099 H 1.8325 -0.2243 H -0.2243 1.8325
0.0000 -0.0001 0.0002
