2
#include<libipv1/ip_lib.h>
4
#include<libciomr/libciomr.h>
5
#include<libint/libint.h>
18
if (Molecule.num_atoms == 0)
21
grad_enuc = block_matrix(Molecule.num_atoms,3);
22
for(i=0; i<Molecule.num_atoms; i++)
23
for(j=0; j<Molecule.num_atoms; j++)
25
rij2 = (Molecule.centers[i].x-Molecule.centers[j].x)*(Molecule.centers[i].x-Molecule.centers[j].x);
26
rij2 += (Molecule.centers[i].y-Molecule.centers[j].y)*(Molecule.centers[i].y-Molecule.centers[j].y);
27
rij2 += (Molecule.centers[i].z-Molecule.centers[j].z)*(Molecule.centers[i].z-Molecule.centers[j].z);
28
tmp = Molecule.centers[i].Z_nuc*Molecule.centers[j].Z_nuc/(rij2*sqrt(rij2));
29
grad_enuc[i][0] -= (Molecule.centers[i].x-Molecule.centers[j].x)*tmp;
30
grad_enuc[i][1] -= (Molecule.centers[i].y-Molecule.centers[j].y)*tmp;
31
grad_enuc[i][2] -= (Molecule.centers[i].z-Molecule.centers[j].z)*tmp;
34
if (UserOptions.print_lvl >= PRINT_OEDERIV)
35
print_atomvec("Nuclear repulsion component of the forces (a.u.)",grad_enuc);
37
add_mat(Grad,grad_enuc,Grad,Molecule.num_atoms,3);
38
free_block(grad_enuc);