1 /*
2  * Copyright 2009-2021 The VOTCA Development Team (http://www.votca.org)
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *     http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  *
16  */
17 
18 #ifndef VOTCA_CSG_CSG_REUPDATE_H
19 #define VOTCA_CSG_CSG_REUPDATE_H
20 
21 // Third party includes
22 #include <boost/program_options.hpp>
23 
24 // VOTCA includes
25 #include <votca/tools/histogramnew.h>
26 #include <votca/tools/property.h>
27 #include <votca/tools/table.h>
28 
29 // Local VOTCA includes
30 #include "votca/csg/csgapplication.h"
31 #include "votca/csg/potentialfunctions/potentialfunction.h"
32 #include "votca/csg/potentialfunctions/potentialfunctioncbspl.h"
33 #include "votca/csg/potentialfunctions/potentialfunctionlj126.h"
34 #include "votca/csg/potentialfunctions/potentialfunctionljg.h"
35 #include "votca/csg/topologyreader.h"
36 
37 using namespace votca::csg;
38 using namespace votca::tools;
39 
40 struct PotentialInfo {
41 
42   PotentialInfo(votca::Index index, bool bonded_, votca::Index vec_pos_,
43                 std::string &param_in_ext_, Property *options,
44                 bool gentable = false);
45 
46   votca::Index potentialIndex;
47   bool bonded;
48   PotentialFunction *ucg;
49   votca::Index vec_pos;
50   std::pair<votca::Index, votca::Index> beadTypes;
51 
52   std::string potentialName;
53   std::string potentialFunction;
54   std::string type1, type2;
55 
56   double rmin, rcut;
57 
58   Property *options_;
59 };
60 
61 class CsgREupdate : public CsgApplication {
62  public:
ProgramName()63   std::string ProgramName() override { return "csg_reupdate"; }
HelpText(std::ostream & out)64   void HelpText(std::ostream &out) override {
65     out << "computes relative entropy update.";
66   }
67 
DoTrajectory()68   bool DoTrajectory() override { return true; }
69 
DoMapping()70   bool DoMapping() override { return false; }
71 
DoThreaded()72   bool DoThreaded() override { return true; }
SynchronizeThreads()73   bool SynchronizeThreads() override { return false; }
74 
NeedsTopology()75   bool NeedsTopology() override { return false; }
76 
77   void Initialize() override;
78   bool EvaluateOptions() override;
79   void BeginEvaluate(Topology *top, Topology *top_atom = nullptr) override;
80   void LoadOptions(const std::string &file);
81 
82   void Run() override;
83 
84   void EndEvaluate() override;
85   std::unique_ptr<CsgApplication::Worker> ForkWorker(void) override;
86   void MergeWorker(Worker *worker) override;
87 
88  private:
89  protected:
90   Property options_;
91   std::vector<Property *> nonbonded_;
92 
93   using PotentialContainer = std::vector<PotentialInfo *>;
94   PotentialContainer potentials_;
95 
96   votca::Index nlamda_;
97   Eigen::VectorXd lamda_;
98   //  HS_ is a symmetric matrix
99   Eigen::MatrixXd HS_;
100   Eigen::VectorXd DS_;
101   Eigen::VectorXd dUFrame_;
102   bool hessian_check_;
103 
104   double UavgAA_;
105   double UavgCG_;
106   double beta_;
107   double relax_;
108   votca::Index nframes_;
109 
110   bool gentable_;
111   bool dosteep_;
112 
113   std::vector<Table *> aardfs_;
114   std::vector<double *> aardfnorms_;
115 
116   // file extension for the inputs/outputs
117   std::string param_in_ext_, param_out_ext_;
118   std::string pot_out_ext_;
119   std::string rdf_ext_;
120 
121   void WriteOutFiles();
122   void EvalBonded(Topology *conf, PotentialInfo *potinfo);
123   void EvalNonbonded(Topology *conf, PotentialInfo *potinfo);
124 
125   // Compute Avg U, dU, and d2U values in reference AA ensemble
126   void AAavgBonded(PotentialInfo *potinfo);
127   void AAavgNonbonded(PotentialInfo *potinfo);
128 
129   // Formulates  HS_ dlamda = -  DS_ system of Lin Eq.
130   void REFormulateLinEq();
131 
132   // Solve  HS_ dlamda = -  DS_ and update  lamda_
133   void REUpdateLamda();
134 };
135 
136 class CsgREupdateWorker : public CsgApplication::Worker {
137  public:
138   ~CsgREupdateWorker() override = default;
139 
140   Property options_;
141   std::vector<Property *> nonbonded_;
142 
143   using PotentialContainer = std::vector<PotentialInfo *>;
144   PotentialContainer potentials_;
145 
146   votca::Index nlamda_;
147   Eigen::VectorXd lamda_;
148   Eigen::MatrixXd HS_;
149   Eigen::VectorXd DS_;
150   Eigen::VectorXd dUFrame_;
151 
152   double UavgCG_;
153   double beta_;
154   votca::Index nframes_;
155 
156   void EvalConfiguration(Topology *conf, Topology *conf_atom) override;
157   void EvalBonded(Topology *conf, PotentialInfo *potinfo);
158   void EvalNonbonded(Topology *conf, PotentialInfo *potinfo);
159 };
160 
161 #endif  // VOTCA_CSG_CSG_REUPDATE_H
162