1 /*-------------------------------------------------------------------
2 Copyright 2013 Ravishankar Sundararaman, Deniz Gunceler
3 
4 This file is part of JDFTx.
5 
6 JDFTx is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10 
11 JDFTx is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 GNU General Public License for more details.
15 
16 You should have received a copy of the GNU General Public License
17 along with JDFTx.  If not, see <http://www.gnu.org/licenses/>.
18 -------------------------------------------------------------------*/
19 
20 #ifndef JDFTX_ELECTRONIC_PCM_H
21 #define JDFTX_ELECTRONIC_PCM_H
22 
23 #include <fluid/FluidSolver.h>
24 #include <core/RadialFunction.h>
25 #include <core/EnergyComponents.h>
26 #include <core/Coulomb.h>
27 
28 //! @addtogroup Solvation
29 //! @{
30 
31 //! Base class for all PCMs
32 class PCM : public FluidSolver
33 {
34 public:
35 	PCM(const Everything& e, const FluidSolverParams& fsp);
36 	virtual ~PCM();
37 
38 	void dumpDensities(const char* filenamePattern) const; //!< dump cavity shape functions
39 	void dumpDebug(const char* filenamePattern) const; //!< generate fluidDebug text file with common info to all PCMs
40 
41 protected:
42 	EnergyComponents Adiel; //!< PCM energy components
43 	ScalarFieldTilde rhoExplicitTilde; //!< Charge density of explicit (electronic) system
44 	ScalarField nCavity, tauCavity, nCavityEx[2]; //!< Cavity determining electron density (or product for SaLSA, or KE density for SG14tauVW, and expanded electron densities for the SGA13 variant)
45 	ScalarFieldArray shape; //!< Electrostatic cavity shape function. Second component, if any, is separate ionic cavity
46 	ScalarField shapeVdw; //!< Separate cavitation/dispersion shape function for the SGA13 variant
47 
printDebug(FILE * fp)48 	virtual void printDebug(FILE* fp) const {} //!< over-ride to get extra PCM-specific output in fluidDebug text file
49 
50 	void updateCavity(); //!< update shape function(s) from nCavity, and energies dependent upon shape alone
51 
52 	//! Propagate A_shape (+ cached Acavity_shape) and accumulate to gradients w.r.t nCavity and rhoExplicitTilde.
53 	//! Set fluid force contributions (for atom-sphere cavities) and stress contributions if non-null.
54 	void propagateCavityGradients(const ScalarFieldArray& A_shape, ScalarField& A_nCavity, ScalarFieldTilde& A_rhoExplicitTilde, IonicGradient* forces, matrix3<>* Adiel_RRT) const;
55 
56 	//! Accumulate extra fluid forces (vdw and full-core forces, when applicable)
57 	void accumExtraForces(IonicGradient* forces, const ScalarFieldTilde& A_nCavityTilde) const;
58 
59 	ScalarFieldTilde getFullCore() const; //!< get full core correction for PCM variants that need them
60 private:
61 	ScalarField Acavity_shape, Acavity_shapeVdw; //!< Cached gradients of cavitation (and dispersion) energies w.r.t shape functions (assumed Acavity does not depend on ionic cavity)
62 	matrix3<> Acavity_RRT; //!< Cached gradients of cavitation (and dispersion) energies w.r.t lattice vectors
63 	double A_nc, A_tension, A_vdwScale, A_eta_wDiel, A_pCavity, A_cavityScale; //!< Cached derivatives w.r.t fit parameters (accessed via dumpDebug() for PCM fits)
64 	double Rex[2]; //!< radii for cavity expansion (SGA13 only)
65 	RadialFunctionG wExpand[2]; //!< weight function for cavity expansion (SGA13 only)
66 	RadialFunctionG wCavity; //!< weight function for nonlocal cavitation energy
67 	std::vector<vector3<>> atposAll; //!< all solute atomic positions (SoftSphere only)
68 	std::vector<vector3<int>> latticeReps; //!< lattice repetitions needed to treat periodic boundaries correctly for SoftSphere
69 	std::vector<double> Rall; //!< all solute atomic radii, with overall scale factor pre-multiplied (SoftSphere only)
70 	std::vector<double> RallIonic; //!< Rall + ion spacing adde, used for separate ionic cavity, if any (SoftSphere only)
71 	ScalarFieldArray zMask; //optional cavity mask function
72 	int nShape; //natural number of shape functions of the solvation model (2 if ionspacing is used to make ionic cavity, else 1)
73 	bool fixedCavityMasked; //!< whether mask has already been applied to fixed cavity
74 protected:
75 	std::vector<RadialFunctionG> Sf; //!< spherically-averaged structure factors for each solvent site
76 	std::vector<int> atomicNumbers; //!< atomic number for each solvent site (for dispersion interactions)
coulomb(const ScalarFieldTilde & rho)77 	static ScalarFieldTilde coulomb(const ScalarFieldTilde& rho) { return (-4*M_PI) * Linv(O(rho)); }
coulombStress(const ScalarFieldTilde & X,const ScalarFieldTilde & Y)78 	static matrix3<> coulombStress(const ScalarFieldTilde& X, const ScalarFieldTilde& Y) { return (-4*M_PI) * LinvStress(X, Y); }
79 	friend struct ChargedDefect;
80 };
81 
82 //! @}
83 #endif // JDFTX_ELECTRONIC_PCM_H
84