1 /*
2 * PCMSolver, an API for the Polarizable Continuum Model
3 * Copyright (C) 2020 Roberto Di Remigio, Luca Frediani and contributors.
4 *
5 * This file is part of PCMSolver.
6 *
7 * PCMSolver is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU Lesser General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * PCMSolver is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public License
18 * along with PCMSolver. If not, see <http://www.gnu.org/licenses/>.
19 *
20 * For information on the complete list of contributors to the
21 * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
22 */
23
24 #include <sys/types.h>
25 #include "Input.hpp"
26
27 #include <algorithm>
28 #include <functional>
29 #include <iostream>
30 #include <map>
31 #include <string>
32 #include <vector>
33
34 #include "Config.hpp"
35
36 #include "utils/getkw/Getkw.h"
37 #include <Eigen/Core>
38
39 #include "PCMInput.h"
40 #include "bi_operators/BIOperatorData.hpp"
41 #include "cavity/CavityData.hpp"
42 #include "green/GreenData.hpp"
43 #include "solver/SolverData.hpp"
44 #include "utils/Factory.hpp"
45 #include "utils/Solvent.hpp"
46 #include "utils/Sphere.hpp"
47
48 namespace pcm {
49
Input(const std::string & filename)50 Input::Input(const std::string & filename) {
51 reader(filename);
52 semanticCheck();
53 }
54
Input(const PCMInput & host_input)55 Input::Input(const PCMInput & host_input) {
56 reader(host_input);
57 semanticCheck();
58 }
59
reader(const std::string & filename)60 void Input::reader(const std::string & filename) {
61 Getkw input_ = Getkw(filename, false, true);
62 units_ = input_.getStr("UNITS");
63 CODATAyear_ = input_.getInt("CODATA");
64 initBohrToAngstrom(bohrToAngstrom, CODATAyear_);
65
66 const Section & mol = input_.getSect("MOLECULE");
67 MEPfromMolecule_ = true;
68 if (mol.isDefined()) {
69 geometry_ = mol.getDblVec("GEOMETRY");
70 MEPfromMolecule_ = mol.getBool("MEP");
71 }
72
73 const Section & cavity = input_.getSect("CAVITY");
74
75 cavityType_ = cavity.getStr("TYPE");
76 area_ = cavity.getDbl("AREA");
77 if (cavityType_ == "RESTART") {
78 cavFilename_ = cavity.getStr("NPZFILE");
79 }
80
81 scaling_ = cavity.getBool("SCALING");
82 radiiSet_ = detail::uppercase(cavity.getStr("RADIISET"));
83 minimalRadius_ = cavity.getDbl("MINRADIUS");
84 mode_ = detail::uppercase(cavity.getStr("MODE"));
85 if (mode_ == "EXPLICIT") {
86 std::vector<double> spheresInput = cavity.getDblVec("SPHERES");
87 int j = 0;
88 int nAtoms = int(spheresInput.size() / 4);
89 for (int i = 0; i < nAtoms; ++i) {
90 Eigen::Vector3d center;
91 center = (Eigen::Vector3d() << spheresInput[j],
92 spheresInput[j + 1],
93 spheresInput[j + 2])
94 .finished();
95 Sphere sph(center, spheresInput[j + 3]);
96 spheres_.push_back(sph);
97 j += 4;
98 }
99 // Initialize molecule from spheres only when molecule section is absent
100 if (!mol.isDefined())
101 molecule_ = Molecule(spheres_);
102 } else if (mode_ == "ATOMS") {
103 atoms_ = cavity.getIntVec("ATOMS");
104 radii_ = cavity.getDblVec("RADII");
105 }
106
107 // Get the contents of the Medium section
108 const Section & medium = input_.getSect("MEDIUM");
109 // Get the name of the solvent
110 std::string name = medium.getStr("SOLVENT");
111 if (name == "EXPLICIT") {
112 hasSolvent_ = false;
113 // Get the probe radius
114 probeRadius_ = medium.getDbl("PROBERADIUS");
115 // Get the contents of the Green<inside> section...
116 const Section & inside = medium.getSect("GREEN<INSIDE>");
117 // ...and initialize the data members
118 greenInsideType_ = inside.getStr("TYPE") + "_" + inside.getStr("DER");
119 epsilonInside_ = inside.getDbl("EPS");
120 // Get the contents of the Green<outside> section...
121 const Section & outside = medium.getSect("GREEN<OUTSIDE>");
122 // ...and initialize the data members
123 greenOutsideType_ = outside.getStr("TYPE") + "_" + outside.getStr("DER");
124 epsilonStaticOutside_ = outside.getDbl("EPS");
125 epsilonDynamicOutside_ = outside.getDbl("EPSDYN");
126 epsilonStatic1_ = outside.getDbl("EPS1");
127 epsilonDynamic1_ = outside.getDbl("EPSDYN1");
128 epsilonStatic2_ = outside.getDbl("EPS2");
129 epsilonDynamic2_ = outside.getDbl("EPSDYN2");
130 center_ = outside.getDbl("CENTER");
131 width_ = outside.getDbl("WIDTH");
132 origin_ = outside.getDblVec("INTERFACEORIGIN");
133 if (outside.getStr("TYPE") == "SPHERICALDIFFUSE") {
134 greenOutsideType_ += "_" + outside.getStr("PROFILE");
135 }
136 maxL_ = outside.getInt("MAXL");
137 } else { // This part must be reviewed!! Some data members are not initialized...
138 // Just initialize the solvent object in this class
139 hasSolvent_ = true;
140 if (solvents().find(name) == solvents().end()) {
141 PCMSOLVER_ERROR("Solvent " + name + " NOT found!");
142 } else {
143 solvent_ = solvents()[name];
144 }
145 probeRadius_ = solvent_.probeRadius * angstromToBohr();
146 // Specification of the solvent by name means isotropic PCM
147 // We have to initialize the Green's functions data here, Solvent class
148 // is an helper class and should not be used in the core classes.
149 greenInsideType_ = "VACUUM_DERIVATIVE";
150 epsilonInside_ = 1.0;
151 greenOutsideType_ = "UNIFORMDIELECTRIC_DERIVATIVE";
152 epsilonStaticOutside_ = solvent_.epsStatic;
153 epsilonDynamicOutside_ = solvent_.epsDynamic;
154 }
155 integratorType_ = medium.getStr("DIAGONALINTEGRATOR");
156 integratorScaling_ = medium.getDbl("DIAGONALSCALING");
157
158 solverType_ = medium.getStr("SOLVERTYPE");
159 correction_ = medium.getDbl("CORRECTION");
160 hermitivitize_ = medium.getBool("MATRIXSYMM");
161 isDynamic_ = medium.getBool("NONEQUILIBRIUM");
162
163 const Section & chgdist = input_.getSect("CHARGEDISTRIBUTION");
164 MEPfromChargeDist_ = false;
165 if (chgdist.isDefined()) {
166 MEPfromChargeDist_ = true;
167 // Set monopoles
168 if (chgdist.getKey<std::vector<double> >("MONOPOLES").isDefined()) {
169 std::vector<double> mono = chgdist.getDblVec("MONOPOLES");
170 int j = 0;
171 int n = int(mono.size() / 4);
172 multipoles_.monopoles = Eigen::VectorXd::Zero(n);
173 multipoles_.monopolesSites = Eigen::Matrix3Xd::Zero(3, n);
174 for (int i = 0; i < n; ++i) {
175 multipoles_.monopolesSites.col(i) =
176 (Eigen::Vector3d() << mono[j], mono[j + 1], mono[j + 2]).finished();
177 multipoles_.monopoles(i) = mono[j + 3];
178 j += 4;
179 }
180 }
181 // Set dipoles
182 if (chgdist.getKey<std::vector<double> >("DIPOLES").isDefined()) {
183 std::vector<double> dipo = chgdist.getDblVec("DIPOLES");
184 int j = 0;
185 int n = int(dipo.size() / 6);
186 multipoles_.dipoles = Eigen::Matrix3Xd::Zero(3, n);
187 multipoles_.dipolesSites = Eigen::Matrix3Xd::Zero(3, n);
188 for (int i = 0; i < n; ++i) {
189 multipoles_.dipolesSites.col(i) =
190 (Eigen::Vector3d() << dipo[j], dipo[j + 1], dipo[j + 2]).finished();
191 multipoles_.dipoles.col(i) =
192 (Eigen::Vector3d() << dipo[j + 3], dipo[j + 4], dipo[j + 5]).finished();
193 j += 6;
194 }
195 }
196 }
197 // Set fluctuating charges
198 isFQ_ = false;
199 const Section & mmfq = input_.getSect("MMFQ");
200 if (mmfq.isDefined()) {
201 isFQ_ = true;
202 isNonPolarizable_ = mmfq.getBool("NONPOLARIZABLE");
203 fragments_.nSitesPerFragment =
204 static_cast<PCMSolverIndex>(mmfq.getInt("SITESPERFRAGMENT"));
205 std::vector<double> sites = mmfq.getDblVec("SITES");
206 int j = 0;
207 int n = int(sites.size() / 5);
208 fragments_.nFragments = int(n / fragments_.nSitesPerFragment);
209 fragments_.sites = Eigen::Matrix3Xd::Zero(3, n);
210 fragments_.chi = Eigen::VectorXd::Zero(n);
211 fragments_.eta = Eigen::VectorXd::Zero(n);
212 for (int i = 0; i < n; ++i) {
213 fragments_.sites.col(i) =
214 (Eigen::Vector3d() << sites[j], sites[j + 1], sites[j + 2]).finished();
215 fragments_.chi(i) = sites[j + 3];
216 fragments_.eta(i) = sites[j + 4];
217 j += 5;
218 }
219 }
220
221 providedBy_ = std::string("API-side");
222 }
223
reader(const PCMInput & host_input)224 void Input::reader(const PCMInput & host_input) {
225 CODATAyear_ = 2010;
226 initBohrToAngstrom(bohrToAngstrom, CODATAyear_);
227
228 cavityType_ = detail::trim_and_upper(host_input.cavity_type);
229 area_ = host_input.area * angstrom2ToBohr2();
230 if (cavityType_ == "RESTART") {
231 cavFilename_ = detail::trim(host_input.restart_name); // No case conversion here!
232 }
233
234 scaling_ = host_input.scaling;
235 radiiSet_ = detail::trim_and_upper(host_input.radii_set);
236 if (radiiSet_ == "UFF") {
237 radiiSetName_ = "UFF";
238 } else if (radiiSet_ == "BONDI") {
239 radiiSetName_ = "Bondi-Mantina";
240 } else {
241 radiiSetName_ = "Allinger's MM3";
242 }
243 minimalRadius_ = host_input.min_radius * angstromToBohr();
244 mode_ = std::string("IMPLICIT");
245
246 std::string name = detail::trim_and_upper(host_input.solvent);
247 if (name.empty() || name == "EXPLICIT") {
248 hasSolvent_ = false;
249 // Get the probe radius
250 probeRadius_ = host_input.probe_radius * angstromToBohr();
251 // Get the contents of the Green<inside> section...
252 // ...and initialize the data members
253 greenInsideType_ =
254 detail::trim_and_upper(host_input.inside_type) + "_DERIVATIVE";
255 epsilonInside_ = 1.0;
256 // Get the contents of the Green<outside> section...
257 // ...and initialize the data members
258 greenOutsideType_ =
259 detail::trim_and_upper(host_input.outside_type) + "_DERIVATIVE";
260 epsilonStaticOutside_ = host_input.outside_epsilon;
261 epsilonDynamicOutside_ = host_input.outside_epsilon;
262 // Initialize interface parameters with bogus values
263 epsilonStatic1_ = 0.0;
264 epsilonDynamic1_ = 0.0;
265 epsilonStatic2_ = 0.0;
266 epsilonDynamic2_ = 0.0;
267 center_ = 0.0;
268 width_ = 0.0;
269 origin_ = std::vector<double>(3, 0.0);
270 maxL_ = 0;
271 } else { // This part must be reviewed!! Some data members are not initialized...
272 // Just initialize the solvent object in this class
273 hasSolvent_ = true;
274 solvent_ = solvents()[name];
275 probeRadius_ = solvent_.probeRadius * angstromToBohr();
276 // Specification of the solvent by name means isotropic PCM
277 // We have to initialize the Green's functions data here, Solvent class
278 // is an helper class and should not be used in the core classes.
279 greenInsideType_ = std::string("VACUUM_DERIVATIVE");
280 epsilonInside_ = 1.0;
281 greenOutsideType_ = std::string("UNIFORMDIELECTRIC_DERIVATIVE");
282 epsilonStaticOutside_ = solvent_.epsStatic;
283 epsilonDynamicOutside_ = solvent_.epsDynamic;
284 }
285
286 integratorType_ = "COLLOCATION";
287 integratorScaling_ = 1.07;
288
289 solverType_ = detail::trim_and_upper(host_input.solver_type);
290 correction_ = host_input.correction;
291 hermitivitize_ = true;
292 isDynamic_ = false;
293 isFQ_ = false;
294
295 providedBy_ = std::string("host-side");
296 }
297
semanticCheck()298 void Input::semanticCheck() {}
299
initMolecule(bool deferred_init)300 void Input::initMolecule(bool deferred_init) {
301 // Gather information necessary to build molecule_
302 // 1. number of atomic centers
303 int nuclei = int(geometry_.size() / 4);
304 // 2. position and charges of atomic centers
305 Eigen::Matrix3Xd centers = Eigen::Matrix3Xd::Zero(3, nuclei);
306 Eigen::VectorXd charges = Eigen::VectorXd::Zero(nuclei);
307 int j = 0;
308 for (int i = 0; i < nuclei; ++i) {
309 centers.col(i) =
310 (Eigen::Vector3d() << geometry_[j], geometry_[j + 1], geometry_[j + 2])
311 .finished();
312 charges(i) = geometry_[j + 3];
313 j += 4;
314 }
315 // 3. list of atoms and list of spheres
316 std::vector<Atom> radiiSet;
317 std::vector<Atom> atoms;
318 atoms.reserve(nuclei);
319 // FIXME Code duplication in function initMolecule in interface/Meddle.cpp
320 tie(radiiSetName_, radiiSet) = utils::bootstrapRadiiSet().create(radiiSet_);
321 for (int i = 0; i < charges.size(); ++i) {
322 int index = int(charges(i)) - 1;
323 atoms.push_back(radiiSet[index]);
324 if (scaling_)
325 atoms[i].radiusScaling = 1.2;
326 }
327 // Based on the creation mode (Implicit or Atoms)
328 // the spheres list might need postprocessing
329 if (mode_ == "IMPLICIT" || mode_ == "ATOMS") {
330 for (int i = 0; i < charges.size(); ++i) {
331 // Convert to Bohr and multiply by scaling factor (alpha)
332 double radius = atoms[i].radius * angstromToBohr() * atoms[i].radiusScaling;
333 spheres_.push_back(Sphere(centers.col(i), radius));
334 }
335 if (mode_ == "ATOMS") {
336 // Loop over the atomsInput array to get which atoms will have a user-given
337 // radius
338 for (size_t i = 0; i < atoms_.size(); ++i) {
339 int index =
340 atoms_[i] - 1; // -1 to go from human readable to machine readable
341 // Put the new Sphere in place of the implicit-generated one
342 spheres_[index] = Sphere(centers.col(index), radii_[i]);
343 }
344 }
345 }
346
347 // 4. masses
348 Eigen::VectorXd masses = Eigen::VectorXd::Zero(nuclei);
349 for (int i = 0; i < masses.size(); ++i) {
350 masses(i) = atoms[i].mass;
351 }
352 // 5. molecular point group
353 // FIXME currently hardcoded to C1
354
355 // OK, now get molecule_
356 molecule_ = Molecule(nuclei, charges, masses, centers, atoms, spheres_);
357 // Check that all atoms have a radius attached
358 std::vector<Atom>::const_iterator res =
359 std::find_if(atoms.begin(), atoms.end(), invalid);
360 if (res != atoms.end() && !deferred_init) {
361 std::cout << molecule_ << std::endl;
362 PCMSOLVER_ERROR("Some atoms do not have a radius attached. Please specify a "
363 "radius for all atoms (see "
364 "http://pcmsolver.readthedocs.org/en/latest/users/input.html)!");
365 }
366 }
367
cavityParams() const368 CavityData Input::cavityParams() const {
369 return CavityData(
370 cavityType_, molecule_, area_, probeRadius_, minimalRadius_, cavFilename_);
371 }
372
insideGreenParams() const373 GreenData Input::insideGreenParams() const {
374 return GreenData(greenInsideType_, epsilonInside_);
375 }
376
outsideStaticGreenParams() const377 GreenData Input::outsideStaticGreenParams() const {
378 GreenData retval(greenOutsideType_, epsilonStaticOutside_);
379 if (not hasSolvent_) {
380 retval.epsilon1 = epsilonStatic1_;
381 retval.epsilon2 = epsilonStatic2_;
382 retval.center = center_;
383 retval.width = width_;
384 retval.origin =
385 (Eigen::Vector3d() << origin_[0], origin_[1], origin_[2]).finished();
386 retval.maxL = maxL_;
387 }
388 return retval;
389 }
390
outsideDynamicGreenParams() const391 GreenData Input::outsideDynamicGreenParams() const {
392 GreenData retval(greenOutsideType_, epsilonDynamicOutside_);
393 if (not hasSolvent_) {
394 retval.epsilon1 = epsilonDynamic1_;
395 retval.epsilon2 = epsilonDynamic2_;
396 retval.center = center_;
397 retval.width = width_;
398 retval.origin =
399 (Eigen::Vector3d() << origin_[0], origin_[1], origin_[2]).finished();
400 retval.maxL = maxL_;
401 }
402 return retval;
403 }
404
solverParams() const405 SolverData Input::solverParams() const {
406 return SolverData(solverType_, correction_, hermitivitize_);
407 }
408
integratorParams() const409 BIOperatorData Input::integratorParams() const {
410 return BIOperatorData(integratorType_, integratorScaling_);
411 }
412
413 namespace detail {
left_trim(std::string s)414 std::string left_trim(std::string s) {
415 s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](int ch) {
416 return !std::isspace(ch);
417 }));
418 return s;
419 }
420
right_trim(std::string s)421 std::string right_trim(std::string s) {
422 s.erase(
423 std::find_if(s.rbegin(), s.rend(), [](int ch) { return !std::isspace(ch); })
424 .base(),
425 s.end());
426 return s;
427 }
428
left_trim(const char * src)429 std::string left_trim(const char * src) {
430 std::string tmp(src);
431 return left_trim(tmp);
432 }
433
right_trim(const char * src)434 std::string right_trim(const char * src) {
435 std::string tmp(src);
436 return right_trim(tmp);
437 }
438
trim(std::string s)439 std::string trim(std::string s) { return right_trim(left_trim(s)); }
440
trim(const char * src)441 std::string trim(const char * src) {
442 std::string tmp(src);
443 return trim(tmp);
444 }
445
uppercase(std::string s)446 std::string uppercase(std::string s) {
447 std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) {
448 return std::toupper(c);
449 });
450 return s;
451 }
452
uppercase(const char * src)453 std::string uppercase(const char * src) {
454 std::string tmp(src);
455 return uppercase(tmp);
456 }
457
trim_and_upper(const char * src)458 std::string trim_and_upper(const char * src) {
459 std::string tmp(src);
460 return uppercase(trim(tmp));
461 }
462
string_to_enum(std::string field)463 PCMInputFields string_to_enum(std::string field) {
464 static std::map<std::string, PCMInputFields> mapper;
465
466 mapper["cavity_type"] = PCMInputFields::cavity_type;
467 mapper["patch_level"] = PCMInputFields::patch_level;
468 mapper["coarsity"] = PCMInputFields::coarsity;
469 mapper["area"] = PCMInputFields::area;
470 mapper["radii_set"] = PCMInputFields::radii_set;
471 mapper["min_distance"] = PCMInputFields::min_distance;
472 mapper["der_order"] = PCMInputFields::der_order;
473 mapper["scaling"] = PCMInputFields::scaling;
474 mapper["restart_name"] = PCMInputFields::restart_name;
475 mapper["min_radius"] = PCMInputFields::min_radius;
476 mapper["solver_type"] = PCMInputFields::solver_type;
477 mapper["correction"] = PCMInputFields::correction;
478 mapper["solvent"] = PCMInputFields::solvent;
479 mapper["probe_radius"] = PCMInputFields::probe_radius;
480 mapper["equation_type"] = PCMInputFields::equation_type;
481 mapper["inside_type"] = PCMInputFields::inside_type;
482 mapper["outside_epsilon"] = PCMInputFields::outside_epsilon;
483 mapper["outside_type"] = PCMInputFields::outside_type;
484
485 return mapper[field];
486 }
487 } // namespace detail
488 } // namespace pcm
489