writing wfks, dens. to restart MD

[gpanchap]

Dear Francois,

Is there a way to write wfks or dens at the end of an MD run (for the final config) to help a smooth restart ? Presently I find that upon restart the energy (and derivatives like stress) is not converged in the first few MD steps.

Thanks,
Ganesh

[fgygi]

The normal procedure to restart a calculation is to save the sample at the end of a run using the "save" command

[qbox] save file.xml

This saves all the information defining the sample on the file "file.xml", including atomic position, cell size and wavefunctions. When restarting, use the "load" command

[qbox] load file.xml

This restores the sample in the same state as it was at the end of the previous run. Note that the charge density is recomputed from the wavefunctions (it is not stored separately in file.xml). In some cases, when using charge mixing between SCF steps, the wavefunctions can be slightly inconsistent with the density (if the SCF iterations are not fully converged). This can lead to a sudden change of the energy upon restarting. This can be avoided by adding an extra run command before saving the sample (with no ionic motion) e.g. "run 0 5 5" (note that 5 may not be sufficient depending on the system studied).

Please post the input file so I can check if there is another problem.

[gpanchap]

I was including a run command at the beginning to converge the energy before proceeding with the MD. I am simulating an insulating solid. Here is the input I used:

#ground state calculation
load pt135NPT.xml.01
set nrowmax 88 #Min (ncpu,z-FFT size)
set ecut 120
set wf_dyn PSD #ppreconditioning
set ecutprec 30 #value for preconditioning
randomize_wf #randomize initial wfs. to avoid saddle points
run 0 500 #initial GS

# start NPT simulation
set thermostat LOWE
set th_temp 300
set th_time 5
set atoms_dyn MD

set stress ON
set ext_stress 0 0 0 0.0 0.0 0.0 #0GPa hyd. str.
set cell_mass 10000000 #total at. mass*1823 (mass of a proton)
set cell_dyn SD
set cell_lock S #length can change but not shape
set dt 41 #1fs=41a.u

print stress
print ext_stress
print cell_mass
print cell_dyn
print cell_lock
print atoms_dyn
print dt
print ref_cell

status
reset_vcm

run 1500 10 # noptimization, nscf, niter within each scf.
save pt135NPT.xml

[fgygi]

(I assume in the following that the file pt135NPT.xml.01 contains the sample saved at the end of a prior simulation) When loading a sample saved in a previous simulation, you should not randomize the wavefunctions (this is only necessary the very first time you compute the electronic ground state to avoid high symmetry spurious solutions). Note that the value of ecut is loaded from the sample file, so it need not be reinitialized. I include below some comments about the input file, in particular regarding the use of the LOWE thermostat:

#ground state calculation
load pt135NPT.xml.01
set nrowmax 88 #Min (ncpu,z-FFT size)
set ecut 120

The command "set ecut 120" is not needed (but does no harm).

set wf_dyn PSD #ppreconditioning
set ecutprec 30 #value for preconditioning
randomize_wf #randomize initial wfs. to avoid saddle points

The command "randomize_wf" should be removed.

run 0 500 #initial GS

This run command can be removed too.

# start NPT simulation
set thermostat LOWE
set th_temp 300
set th_time 5

Note that this value of th_time is extremely small and will result in a very high collision frequency. This may significantly affect the dynamics. It would be advisable to test whether larger values (up to several thousand) lead to similar results.

set atoms_dyn MD

set stress ON
set ext_stress 0 0 0 0.0 0.0 0.0 #0GPa hyd. str.

This command is not necessary. The external stress is zero by default.

set cell_mass 10000000 #total at. mass*1823 (mass of a proton)
set cell_dyn SD
set cell_lock S #length can change but not shape
set dt 41 #1fs=41a.u

print stress
print ext_stress
print cell_mass
print cell_dyn
print cell_lock
print atoms_dyn
print dt
print ref_cell

status
reset_vcm

run 1500 10 # noptimization, nscf, niter within each scf.
save pt135NPT.xml