Developers:Forces
From OctopusWiki
Contents |
Forces from the gradient of the wavefunctions
Formulae
Comparisons
Here we compare the forces calculated using either the derivative of the potential or the derivative of the wavefunctions.
Converge of the forces with respect to the grid spacing
Na2 (local pseudopotential)
N2 (non-local pseudopotential)
Eggbox effect
The change in the magnitude of the force, when the position of the molecule is changed inside the box (keeping the geometry of the molecule invariant).
Spacing is 0.5 [b] for all cases.
Na2 (local pseudopotential)
N2 (non-local pseudopotential)
The next graphic is for N2 and shows the difference between the magnitude of the force over the atoms (
). Ideally this should be zero, but is not because of the grid.
Input files
This is the input file used for the N tests (a bond length of 1.5 [b] was used to measure the eggbox effect)
SystemName = "N2" CalculationMode = gs bond_length = 2.0744 %Coordinates "N" | -bond_length/2 + pos | 0.0 | 0.0 | no "N" | bond_length/2 + pos | 0.0 | 0.0 | no % BoxShape = sphere radius = 14 ConvRelDens = 1e-9
and this is the bash script that calls the input file, it contains the definition of the missing variables
inpfile=inp
for OCT_spacing in `seq -w 0.8 -0.05 0.19`
do
for OCT_forces in derivate_wavefunctions derivate_potential
do
for OCT_pos in 0.1 #`seq -w 0 0.05 0.51`
do
export OCT_PARSE_ENV=1
export OCT_forces
export OCT_pos
export OCT_spacing
dir=run,$OCT_forces,$OCT_spacing,$OCT_pos,
mkdir $dir
cd $dir
cp ../$inpfile inp
octopus
cd ..
done
done
done
Performance
Unfortunately, maintaining all other parameters equal, the calculation of the forces by derivating the wavefunctions is more expensive, the reason being that one needs to calculate the gradient of each Kohn-Sham wavefunction. Note that if there are no non-local potential parts, then one only needs to calculate the gradient of the density. However, we will typically use non-local pseudopotentials.
This effect is negligible if one just does a ground state calculation, and is only interested in looking at the forces -- or even if we do geometry optimization, since the time it takes to do the SCF cycle will be much larger than the time it takes to do the final force calculation. Therefore, in those cases it is obvious that the new scheme (wavefunctions derivations) is better. Especially if we finally put a series effort into geometry minimization (about time) this could make a difference.
However, when doing Molecular Dynamics the performance difference could be relevant. Here I paste some quick profiling results.
The input file used is:
ProfilingMode = yes SystemName = "N2" FromScratch = yes #Forces = derivate_potential Forces = derivate_wavefunctions CalculationMode = td bond_length = 2.0744 pos = 0.0 %Coordinates "N" | -bond_length/2 + pos | 0.0 | 0.0 | no "N" | bond_length/2 + pos | 0.0 | 0.0 | no % BoxShape = sphere radius = 14 ConvRelDens = 1e-9 T = 0.020 dt = 0.001 TDMaximumIter = T/dt TDEvolutionMethod = aetrs TDExponentialMethod = taylor TDExpOrder = 4 TDTimeStep = dt AbsorbingBoundaries = no MoveIons = vel_verlet RestartFileFormat = restart_binary
The profiling file obtained after running this with Forces=derivate_potential is:
TAG NUMBER OF CALLS TOTAL TIME TIME PER CALL ======================================================================= COMPLETE_DATASET 1 69.322186 69.3221893 MF_INTEGRATE 256 0.170462 0.0006659 MF_DOTP 153 0.177901 0.0011628 MF_NRM2 2 0.000681 0.0003405 NL_OPERATOR 905 17.221615 0.0190294 HPSI 905 30.021140 0.0331725 KINETIC 905 23.829708 0.0263312 VLPSI 905 4.880140 0.0053924 VNLPSI 905 0.090671 0.0001002 POISSON_SOLVE 23 16.155236 0.7024015 XC 23 2.127298 0.0924912 XC_LOCAL 23 1.964610 0.0854178 TIME_STEP 20 61.513029 3.0756514 DISK_ACCESS 24 2.759057 0.1149607 FORCES 21 1.897243 0.0903449 EPOT_GENERATE 41 6.234250 0.1520549 ELEC_PROPAGATOR 20 31.761460 1.5880730
whereas if we use Forces=derivate_wavefunctions:
TAG NUMBER OF CALLS TOTAL TIME TIME PER CALL ======================================================================= COMPLETE_DATASET 1 74.167379 74.1673813 MF_INTEGRATE 256 0.195562 0.0007639 MF_DOTP 153 0.220761 0.0014429 MF_NRM2 2 0.000739 0.0003695 NL_OPERATOR 1220 19.935452 0.0163405 HPSI 905 30.182227 0.0333505 KINETIC 905 24.080008 0.0266077 VLPSI 905 4.733472 0.0052304 VNLPSI 905 0.092900 0.0001027 POISSON_SOLVE 23 17.153588 0.7458081 XC 23 2.231594 0.0970258 XC_LOCAL 23 2.047600 0.0890261 TIME_STEP 20 66.375619 3.3187809 DISK_ACCESS 24 2.764224 0.1151760 FORCES 21 5.975013 0.2845244 EPOT_GENERATE 41 5.645966 0.1377065 ELEC_PROPAGATOR 20 32.171846 1.6085923
Note that the time it takes to calculate the forces is multiplied roughly by three. increasing the time of the td simulation (TIME_STEP) from 61.5 to 66.3 (before the improvement of the code done by Xavier, the results were much worse: the time it took to calculate the forces was multiplied roughly by eight, increasing the time of the td simulation (TIME_STEP) from 58.3 to 70.
Overall, an increase in a typical MD simulation with moving ions of 5%. This has to be weighted agains a couple of advantages:
- The code gets much cleaner and easier to maintain.
- In principle, the convergence of the forces with respect to grid spacing is much better with the new scheme, which should allow for a smaller grid
spacing, and therefore a gain in speed probably larger than that 5%.
