1/*
2 *                This source code is part of
3 *
4 *                     E  R  K  A  L  E
5 *                             -
6 *                       HF/DFT from Hel
7 *
8 * Written by Susi Lehtola, 2010-2013
9 * Copyright (c) 2010-2013, Susi Lehtola
10 *
11 * This program is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU General Public License
13 * as published by the Free Software Foundation; either version 2
14 * of the License, or (at your option) any later version.
15 */
16
17/* Restricted case */
18#if defined(RESTRICTED) && defined(DFT)
19arma::vec SCF::force_RDFT(rscf_t & sol, const std::vector<double> & occs, const dft_t dft, DFTGrid & grid, DFTGrid & nlgrid, double tol)
20
21#elif defined(RESTRICTED) && defined(HF)
22arma::vec SCF::force_RHF(rscf_t & sol, const std::vector<double> & occs, double tol)
23
24#elif defined(UNRESTRICTED) && defined(DFT)
25arma::vec SCF::force_UDFT(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, const dft_t dft, DFTGrid & grid, DFTGrid & nlgrid, double tol)
26
27#elif defined(UNRESTRICTED) && defined(HF)
28arma::vec SCF::force_UHF(uscf_t & sol, const std::vector<double> & occa, const std::vector<double> & occb, double tol)
29
30#elif defined(UNRESTRICTED) && defined(_ROHF)
31arma::vec SCF::force_ROHF(uscf_t & sol, int Nel_alpha, int Nel_beta, double tol)
32#endif
33{
34
35  arma::mat W;
36#ifdef RESTRICTED
37  W=form_density(sol.E,sol.C,occs);
38#else
39  W=form_density(sol.Ea,sol.Ca,occa)+form_density(sol.Eb,sol.Cb,occb);
40#endif
41
42  // Compute force
43  arma::vec fpul_kin=basisp->kinetic_pulay(sol.P);
44  //    interpret_force(fpul_kin).print("Kinetic Pulay");
45
46  arma::vec fpul_nuc=basisp->nuclear_pulay(sol.P);
47  //    interpret_force(fpul_nuc).print("Nuclear Pulay");
48
49  arma::vec fnuc=basisp->nuclear_der(sol.P);
50  //    interpret_force(fnuc).print("Hellman-Feynman");
51
52  arma::vec forth=basisp->overlap_der(W);
53  //    interpret_force(forth).print("Orthonormality");
54
55  arma::vec frep=basisp->nuclear_force();
56  //    interpret_force(frep).print("Nuclear repulsion");
57
58#ifdef DFT
59  // Get range separation info
60  double omega, kfull, kshort;
61  range_separation(dft.x_func,omega,kfull,kshort);
62#else
63  double kfull=1.0;
64#endif
65
66  // Exact exchange and Coulomb
67  arma::vec fx_full;
68  fx_full.zeros(fnuc.n_elem);
69
70#ifdef DFT
71  // Short-range exchange
72  arma::vec fx_short;
73  fx_short.zeros(fnuc.n_elem);
74
75  if(kfull==0.0) {
76
77    if(densityfit) {
78      if(kshort!=0.0 || kfull!=0.0)
79	throw std::runtime_error("Forces not implemented for density fitting of exact exchange.\n");
80
81      fx_full=dfit.forceJ(sol.P);
82    } else {
83      if(!direct)
84	scr.fill(basisp,intthr,verbose);
85
86      fx_full=scr.forceJ(sol.P,tol);
87    }
88  } else {
89#endif
90
91    if(!direct)
92      scr.fill(basisp,intthr,verbose);
93
94#ifdef RESTRICTED
95    fx_full=scr.forceJK(sol.P,tol,kfull);
96#else
97    fx_full=scr.forceJK(sol.Pa,sol.Pb,tol,kfull);
98#endif
99
100#ifdef DFT
101  }
102
103  if(omega != 0.0) {
104    // Set range separation and fill (low-cost compared to forces anyway)
105    scr_rs.set_range_separation(omega,0.0,1.0);
106    scr_rs.fill(basisp,intthr,verbose);
107
108    // Exchange forces
109#ifdef RESTRICTED
110    fx_short=scr_rs.forceK(sol.P,tol,kshort);
111#else
112    fx_short=scr_rs.forceK(sol.Pa,sol.Pb,tol,kshort);
113#endif
114  }
115
116  // Get the DFT contribution
117  arma::vec fxc;
118  fxc.zeros(fnuc.n_elem);
119  if(dft.x_func>0 || dft.c_func>0) {
120#ifdef RESTRICTED
121    fxc=grid.eval_force(dft.x_func,dft.c_func,sol.P);
122#else
123    fxc=grid.eval_force(dft.x_func,dft.c_func,sol.Pa,sol.Pb);
124#endif
125  }
126
127  // Non-local contribution to force
128  if(dft.nl) {
129    arma::vec vv10f(grid.eval_VV10_force(nlgrid,dft.vv10_b,dft.vv10_C,sol.P));
130    //interpret_force(vv10f).t().print("VV10 force");
131    fxc+=vv10f;
132  }
133
134#endif // DFT
135
136  arma::vec ftot=fpul_kin+fpul_nuc+fnuc+forth+frep;
137  arma::vec ffull=ftot+fx_full;
138#ifdef DFT
139  // DFT contribution
140  ffull+=fxc+fx_short;
141#endif
142
143  if(lincalc) {
144    // Get rid of numerical noise for linear molecule restriction
145    for(size_t inuc=0;inuc<ffull.n_elem/3;inuc++) {
146      ffull(3*inuc)=0.0;
147      ffull(3*inuc+1)=0.0;
148    }
149  }
150
151  //  interpret_force(fx_full).print("Coulomb + exact exchange");
152  //  interpret_force(ftot).print("Total (w.o. 2-electron)");
153  // interpret_force(ffull).t().print("Total force");
154
155  return ffull;
156}
157