# Developers Manual:Forces

## Forces from the gradient of the wavefunctions

### Formulae $f^A_i =-\frac{\partial E_{loc}}{\partial R_i^A} = -2\Re \int\frac{\partial n(r)}{\partial r_i}v_{loc}\left(|r-R^A|\right)dr$

## 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 ( $|f_{atom\ 1}| - |f_{atom\ 2}|$). 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

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

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%.