Surface forces

[moritzgubler]

This is my implementation of the idea I had of calculating forces with surface integrals.

It works quite well and the forces are more accurate. Also, mrchem struggles to compute forces when the world precision is 1e-7 or smaller. Then it requires a lot of memory (more than 32 GB for a H2 molecule). This is not the case with my approach . Geometry optimizations seem to be reliable if the stopping criterion is chosen 10 * world_prec.

I have already put all the parameters into the parser, so it is quite simple to test and run my code (look for the section "- name: Forces" in the template.yaml input parser file).

It works with lda and gga functionals and for both closed and open shell systems.
At the moment, there is one thing that might be improved in the future:

• The exact exchange stress is not computed. Therefore forces can not be computed for exact exchange. It shouldn't be too complicated to implement, but in a quick search I did not find a reference that derives it.

[moritzgubler]

Hi @robertodr @ilfreddy
I implemented now everything I wanted to, cleaned up to code and added a bit of documentation. The new method to calculate the forces now works for LDA and GGA functionals for both restricted and unrestricted calculations. The pulll request is now ready for review.
Best, Moritz

[moritzgubler]

I previously tested GGA calculations only for H2 where the forces where the forces were correct. In more complicated systems this is not the case. (LDA is always correct). I will look into it and let you know once I have solved the problem.

[moritzgubler]

The GGA bug is fixed now, I forgot to add the divergence terms in the xc potential for the GGA case. They have been added now.
MPI does not work anymore, it segfaults in the calculation of the kinetic stress. I suspect, that I try to access orbitals in an orbitalvector that are stored on another rank. @gitpeterwind Do you know what could be the problem?

[moritzgubler]

@gitpeterwind These are the lines that cause the segfault

[gitpeterwind]

In an OrbitalVector, not all orbitals are directly available to all MPI processes. You should do something like:

     if (mrcpp::mpi::my_orb(phi_i)) {
           pos[0] = gridPos(i, 0);
            pos[1] = gridPos(i, 1);
            pos[2] = gridPos(i, 2);
            n1 = nablaPhi[iOrb][0].real().evalf(pos);
            n2 = nablaPhi[iOrb][1].real().evalf(pos);
            n3 = nablaPhi[iOrb][2].real().evalf(pos);
            stress[i](0, 0) -= occ * n1 * n1;
            stress[i](1, 1) -= occ * n2 * n2;
            stress[i](2, 2) -= occ * n3 * n3;
            stress[i](0, 1) -= occ * n1 * n2;
            stress[i](0, 2) -= occ * n1 * n3;
            stress[i](1, 2) -= occ * n2 * n3;
           }
         }
    }

And then use      mrcpp::mpi::allreduce_vector to collect all the results. It is a bit difficult with the way you have defined the stress vector. But if you cannot easily define it as a DoubleVector, or DoubleMatrix , I can try to give a detailed way to do it (one way is simply to copy all the values into a DoubleVector, the do the allreduce operation, and copy back). 

[moritzgubler]

That makes sense, thanks for your quick answer. I will use arrays for the stress tensor.

[moritzgubler]

@gitpeterwind I think it might be a bit more complicated, I need

            nablaPhi[iOrb][0]
            nablaPhi[iOrb][1]
            nablaPhi[iOrb][2]

all on the same rank. How can I make sure that this is the case?

[gitpeterwind]

I think it is correct: the last [0],[1],[2] picks out one OrbitalVector, and iOrb only determines the rank. Or did I misunderstood?

[moritzgubler]

No, nablaPhi[0][0] contains the x component of the gradient of the zeroth orbital. nablaPhi[0] is an orbital vector that contains the x, y and z component of the gradient of orbital 0.

[moritzgubler]

Adding the following print to the loop

    for (int i = 0; i < Phi.size(); i++) {
        nablaPhi.push_back(nabla(Phi[i]));
        nablaPhi[i].distribute();
        if (mrcpp::wrk_rank == 0){
            std::cout << "gradient of orbital " << i << " stored on rank: x: " << nablaPhi[i][0].getRank() << " y: " << nablaPhi[i][1].getRank() << " z: " << nablaPhi[i][2].getRank() << std::endl;
        }
    }

results in:

gradient of orbital 0 stored on rank: x: 0 y: 1 z: 2
gradient of orbital 1 stored on rank: x: 0 y: 1 z: 2
gradient of orbital 2 stored on rank: x: 0 y: 1 z: 2
gradient of orbital 3 stored on rank: x: 0 y: 1 z: 2

but I need all components of the gradient on the same rank

[gitpeterwind]

The OrbitalVector class expects some properties, like that the size is the "usual" number of orbitals, and that each MPI process only access its own indices. If a vector has only have 3 orbitals, you should not define it as an OrbitalVector (and "distribute them is meaningless). Can you simply define std::vector<std::vector<Orbital>> nablaPhi(3); ? or do you need some OrbitalVector specific properties? Still Phi[i] will only be defined on the MPI process with the right rank.
Anyway, if you use the if (mrcpp::mpi::my_orb(phi_i)) { , do you still get the seg fault?

[moritzgubler]

This is a good idea. I thought that this wouldn't work because nabla(Phi) returns an orbitalvector. I tried your suggestion and it compiled. I will try to solve the problem like that, thanks for your help

[gitpeterwind]

There was a bug in the master code. A PR has been sent

[moritzgubler]

Here I create the vector that contains the gradient of all orbitals

[gitpeterwind]

I think it is necessary (at least it won't harm) to add:

std::vector<mrchem::OrbitalVector> nablaPhi;
    mrchem::OrbitalVector hessRho = hess(rho);
    for (int i = 0; i < Phi.size(); i++) {
        nablaPhi.push_back(nabla(Phi[i]));
    }
    nablaPhi.distribute();

( distribute will tell each element in the vector which rank it has)

[moritzgubler]

ok, i will add this

[moritzgubler]

I can recreate the error with this input:

world_prec = 1.0e-4
world_unit = bohr

WaveFunction {
  method = blyp
  restricted = false
}

Properties {
  geometric_derivative = true
}
Forces {
    method = "surface_integrals"
    #method = 'hellmann_feynman'
    surface_integral_precision = medium
}

Molecule {
$coords
H 0.86387418 -0.02789966 0.13501906
Li -0.44463806 -0.04858572 -0.04959119
$end
}

I start the simulation with:

mrchem --dryrun hli.inp
mpirun -np 4 mrchem.x hli.json

[ilfreddy]

@moritzgubler would you mind adding also some testing to your code? According to the code coverage report, your patch is only marginally covered. Tests are essential to make sure the code does not get broken by others later on.

[ilfreddy]

Could you add some class documentation?

[moritzgubler]

@ilfreddy I added some documentation

[ilfreddy]

Just a general comment from my side: could you add some extra testing:

1. There is no unit-testing as far as I can see
2. There is only one regression test

Some classes could also benefit from some additional documentation.

As for the details of the code, I guess it is best if Stig or Peter comment on those.

[moritzgubler]

I added a test, is there a way to look at the codecov report?

[ilfreddy]

This looks good to me now. Just a couple of minor questions.

@moritzgubler

1. could you update the description at the top of the PR to reflect the actual state of the PR (it still says GGA is not included).
2. a question: how about other functionals (like meta-GGA), Would they require additional implementation
3. Is HF exchange working now?

@stigrj this branch includes a merge a few commits ago. Would you let it pass or should it be changed to a rebase?

[moritzgubler]

@ilfreddy

1. lda, gga, closed and open shell are working, the description is updated.
2. the meta-GGA contribution to the stress would have to be implemented. This should be quite easy. I thought only l and gga functionals are implemented in mrchem.
3. No. This would require the exchange stress density. This paper might help: https://pubs.acs.org/doi/10.1021/acs.jctc.0c01194?ref=pdf

[ilfreddy]

@ilfreddy

1. lda, gga, closed and open shell are working, the description is updated.
2. the meta-GGA contribution to the stress would have to be implemented. This should be quite easy. I thought only l and gga functionals are implemented in mrchem.
3. No. This would require the exchange stress density. This paper might help: https://pubs.acs.org/doi/10.1021/acs.jctc.0c01194?ref=pdf

Thanks for the reference! Indeed we only have GGA at present, but we have some plans to make a libXC interface and thus include also metaGGA functionals. But I guess at this point it might be more relevant to focus on HF exchange.