1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman, Kendra Letchworth Weaver, 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 
21 #ifndef JDFTX_ELECTRONIC_NONLINEARPCM_H
22 #define JDFTX_ELECTRONIC_NONLINEARPCM_H
23 
24 #include <fluid/PCM.h>
25 #include <core/VectorField.h>
26 #include <core/Minimize.h>
27 #include <core/Pulay.h>
28 
29 //! @addtogroup Solvation
30 //! @{
31 //! @file NonlinearPCM.h NonlinearPCM and helper classes
32 
33 namespace NonlinearPCMeval { struct Screening; struct Dielectric; } //Forward declaration of helper classes
34 
35 typedef ScalarFieldMultiplet<ScalarFieldData,5> ScalarFieldMuEps; //!< ion chemical potentials and effective local electric field
36 
37 //! Nonlinear solvation models: shared electrostatic part implementation
38 class NonlinearPCM : public PCM, public Minimizable<ScalarFieldMuEps>, public Pulay<ScalarFieldTilde>
39 {
40 public:
41 	ScalarFieldMuEps state; //!< State of the solver = ion chemical potentials and effective local electric field
42 
43 	//! See createFluidSolver()
44 	NonlinearPCM(const Everything& e, const FluidSolverParams& params);
45     virtual ~NonlinearPCM();
46 
prefersGummel()47 	bool prefersGummel() const { return true; }
48 
49 	void loadState(const char* filename); //!< Load state from file
50 	void saveState(const char* filename) const; //!< Save state to file
51 	void dumpDensities(const char* filenamePattern) const;
52 	void minimizeFluid(); //!< Converge using nonlinear conjugate gradients
53 
54 	//! Compute gradient and free energy (with optional outputs)
55 	double operator()(const ScalarFieldMuEps& state, ScalarFieldMuEps& Adiel_state,
56 		ScalarFieldTilde* Adiel_rhoExplicitTilde=0, ScalarFieldTilde* Adiel_nCavityTilde=0, IonicGradient* forces=0, matrix3<>* Adiel_RRT=0) const;
57 
58 	// Interface for Minimizable:
59 	void step(const ScalarFieldMuEps& dir, double alpha);
60 	double compute(ScalarFieldMuEps* grad, ScalarFieldMuEps* Kgrad);
61 
62 protected:
63 	void set_internal(const ScalarFieldTilde& rhoExplicitTilde, const ScalarFieldTilde& nCavityTilde);
64 	double get_Adiel_and_grad_internal(ScalarFieldTilde& Adiel_rhoExplicitTilde, ScalarFieldTilde& Adiel_nCavityTilde, IonicGradient* extraForces, matrix3<>* Adiel_RRT) const;
65 
66 private:
67 	double pMol, ionNbulk, ionZ;
68 	NonlinearPCMeval::Screening* screeningEval; //!< Internal helper class for Screening from PCM_internal
69 	NonlinearPCMeval::Dielectric* dielectricEval; //!< Internal helper class for Dielectric from PCM_internal
70 	RadialFunctionG preconditioner; //!< preconditioner for minimizer version
71 	std::shared_ptr<RealKernel> metric; //!< Pulay metric for SCF version
72 	RadialFunctionG gLookup, xLookup; //!< lookup tables for transcendental solutions involved in the dielectric and ionic SCF method
73 	std::shared_ptr<class LinearPCM> linearPCM;
74 
75 protected:
76 	//Interface for Pulay<ScalarFieldTilde>
77 	double cycle(double dEprev, std::vector<double>& extraValues);
axpy(double alpha,const ScalarFieldTilde & X,ScalarFieldTilde & Y)78 	void axpy(double alpha, const ScalarFieldTilde& X, ScalarFieldTilde& Y) const { ::axpy(alpha, X, Y); }
dot(const ScalarFieldTilde & X,const ScalarFieldTilde & Y)79 	double dot(const ScalarFieldTilde& X, const ScalarFieldTilde& Y) const { return ::dot(X, Y); }
variableSize()80 	size_t variableSize() const { return gInfo.nG * sizeof(complex); }
81 	void readVariable(ScalarFieldTilde& X, FILE* fp) const;
82 	void writeVariable(const ScalarFieldTilde& X, FILE* fp) const;
83 	ScalarFieldTilde getVariable() const;
84 	void setVariable(const ScalarFieldTilde&);
85 	ScalarFieldTilde precondition(const ScalarFieldTilde&) const;
86 	ScalarFieldTilde applyMetric(const ScalarFieldTilde&) const;
87 private:
88 	void phiToState(bool setState); //!< update state if setState=true and epsilon/kappaSq in linearPCM if setState=false from the current phi
89 };
90 
91 //! @}
92 #endif // JDFTX_ELECTRONIC_NONLINEARPCM_H
93