Skip to content

Commit 7bc4631

Browse files
robertodrrobertodr
committed
Add parsing of classical multipolar charge distributions from input
1 parent 770243f commit 7bc4631

10 files changed

Lines changed: 354 additions & 22 deletions

File tree

‎src/interface/Input.cpp‎

Lines changed: 21 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -166,12 +166,12 @@ void Input::reader(const std::string & filename) {
166166
std::vector<double> mono = chgdist.getDblVec("MONOPOLES");
167167
int j = 0;
168168
int n = int(mono.size() / 4);
169-
fragments_.monopoles = Eigen::VectorXd::Zero(n);
170-
fragments_.monopolesSites = Eigen::Matrix3Xd::Zero(3, n);
169+
multipoles_.monopoles = Eigen::VectorXd::Zero(n);
170+
multipoles_.monopolesSites = Eigen::Matrix3Xd::Zero(3, n);
171171
for (int i = 0; i < n; ++i) {
172-
fragments_.monopolesSites.col(i) =
172+
multipoles_.monopolesSites.col(i) =
173173
(Eigen::Vector3d() << mono[j], mono[j + 1], mono[j + 2]).finished();
174-
fragments_.monopoles(i) = mono[j + 3];
174+
multipoles_.monopoles(i) = mono[j + 3];
175175
j += 4;
176176
}
177177
}
@@ -180,31 +180,34 @@ void Input::reader(const std::string & filename) {
180180
std::vector<double> dipo = chgdist.getDblVec("DIPOLES");
181181
int j = 0;
182182
int n = int(dipo.size() / 6);
183-
fragments_.dipoles = Eigen::Matrix3Xd::Zero(3, n);
184-
fragments_.dipolesSites = Eigen::Matrix3Xd::Zero(3, n);
183+
multipoles_.dipoles = Eigen::Matrix3Xd::Zero(3, n);
184+
multipoles_.dipolesSites = Eigen::Matrix3Xd::Zero(3, n);
185185
for (int i = 0; i < n; ++i) {
186-
fragments_.dipolesSites.col(i) =
186+
multipoles_.dipolesSites.col(i) =
187187
(Eigen::Vector3d() << dipo[j], dipo[j + 1], dipo[j + 2]).finished();
188-
fragments_.dipoles.col(i) =
188+
multipoles_.dipoles.col(i) =
189189
(Eigen::Vector3d() << dipo[j + 3], dipo[j + 4], dipo[j + 5]).finished();
190190
j += 6;
191191
}
192192
}
193193
// Set fluctuating charges
194194
isFQ_ = false;
195-
if (chgdist.getKey<std::vector<double> >("FQ").isDefined()) {
195+
const Section & mmfq = input_.getSect("MMFQ");
196+
if (mmfq.isDefined()) {
196197
isFQ_ = true;
197-
std::vector<double> fq = chgdist.getDblVec("FQ");
198+
fragments_.nSitesPerFragment =
199+
static_cast<PCMSolverIndex>(mmfq.getInt("SITESPERFRAGMENT"));
200+
std::vector<double> sites = mmfq.getDblVec("SITES");
198201
int j = 0;
199-
int n = int(fq.size() / 5);
200-
fragments_.FQSites = Eigen::Matrix3Xd::Zero(3, n);
201-
fragments_.FQChi = Eigen::VectorXd::Zero(n);
202-
fragments_.FQEta = Eigen::VectorXd::Zero(n);
202+
int n = int(sites.size() / 5);
203+
fragments_.sites = Eigen::Matrix3Xd::Zero(3, n);
204+
fragments_.chi = Eigen::VectorXd::Zero(n);
205+
fragments_.eta = Eigen::VectorXd::Zero(n);
203206
for (int i = 0; i < n; ++i) {
204-
fragments_.FQSites.col(i) =
205-
(Eigen::Vector3d() << fq[j], fq[j + 1], fq[j + 2]).finished();
206-
fragments_.FQChi(i) = fq[j + 3];
207-
fragments_.FQEta(i) = fq[j + 4];
207+
fragments_.sites.col(i) =
208+
(Eigen::Vector3d() << sites[j], sites[j + 1], sites[j + 2]).finished();
209+
fragments_.chi(i) = sites[j + 3];
210+
fragments_.eta(i) = sites[j + 4];
208211
j += 5;
209212
}
210213
}

‎src/interface/Input.hpp‎

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -46,12 +46,14 @@ struct SolverData;
4646
} // namespace pcm
4747

4848
#include "utils/ChargeDistribution.hpp"
49+
#include "utils/MMFQ.hpp"
4950
#include "utils/Molecule.hpp"
5051
#include "utils/Solvent.hpp"
5152
#include "utils/Sphere.hpp"
5253

5354
namespace pcm {
5455
using utils::ChargeDistribution;
56+
using utils::MMFQ;
5557
using utils::Solvent;
5658
using utils::Sphere;
5759

@@ -131,7 +133,8 @@ class PCMSolver_EXPORT Input {
131133
BIOperatorData integratorParams() const;
132134
/// @}
133135

134-
ChargeDistribution fragments() const { return fragments_; }
136+
ChargeDistribution multipoles() const { return multipoles_; }
137+
MMFQ fragments() const { return fragments_; }
135138
bool MEPfromMolecule() { return MEPfromMolecule_; }
136139

137140
/// Operators
@@ -246,8 +249,10 @@ class PCMSolver_EXPORT Input {
246249
bool isFQ_;
247250
/// Whether to calculate the MEP from the molecular geometry
248251
bool MEPfromMolecule_;
249-
/// Classical charge distribution
250-
ChargeDistribution fragments_;
252+
/// Classical charge distribution of point multipoles
253+
ChargeDistribution multipoles_;
254+
/// Classical fluctuating charges MM force field
255+
MMFQ fragments_;
251256
/// Who performed the syntactic input parsing
252257
std::string providedBy_;
253258
};

‎src/mmfq/CMakeLists.txt‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
# List of headers
2+
list(APPEND headers_list FQOhno.hpp)
3+
4+
# List of sources
5+
list(APPEND sources_list FQOhno.cpp)
6+
7+
add_library(mmfq OBJECT ${sources_list} ${headers_list})
8+
set_target_properties(mmfq PROPERTIES POSITION_INDEPENDENT_CODE 1 )
9+
add_dependencies(mmfq generate-config-hpp)
10+
target_compile_options(mmfq PRIVATE "$<$<CONFIG:DEBUG>:${EXDIAG_CXX_FLAGS}>")
11+
# Sets install directory for all the headers in the list
12+
foreach(_header ${headers_list})
13+
install(FILES ${_header} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}/mmfq)
14+
endforeach()

‎src/mmfq/FQOhno.cpp‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,83 @@
1+
/**
2+
* PCMSolver, an API for the Polarizable Continuum Model
3+
* Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.
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 "FQOhno.hpp"
25+
26+
#include <iostream>
27+
28+
#include "Config.hpp"
29+
30+
#include <Eigen/Cholesky>
31+
#include <Eigen/Core>
32+
33+
#include "utils/MMFQ.hpp"
34+
35+
namespace pcm {
36+
namespace mmfq {
37+
using utils::MMFQ;
38+
39+
FQOhno::FQOhno(const MMFQ & ff) : mmfq_(ff) { buildSystemMatrix_impl(); }
40+
41+
void FQOhno::buildSystemMatrix_impl() {
42+
PCMSolverIndex nFragments = mmfq_.nFragments;
43+
PCMSolverIndex nSitesPerFragment = mmfq_.nSitesPerFragment;
44+
PCMSolverIndex nSites = nFragments * nSitesPerFragment;
45+
PCMSolverIndex fq_dim = nSites + nFragments;
46+
TIMER_ON("Computing Dlambda");
47+
Dlambda_ = Eigen::MatrixXd::Zero(fq_dim, fq_dim);
48+
// Fill J block
49+
for (PCMSolverIndex i = 0; i < nSites; ++i) {
50+
for (PCMSolverIndex j = 0; j < nSites; ++j) {
51+
if (i == j) {
52+
Dlambda_(i, i) = mmfq_.eta(i);
53+
} else {
54+
double eta_ij = 0.5 * (mmfq_.eta(i) + mmfq_.eta(j));
55+
double dist = (mmfq_.sites.col(i) - mmfq_.sites.col(j)).norm();
56+
Dlambda_(i, j) = eta_ij / std::sqrt(1.0 + std::pow(eta_ij * dist, 2));
57+
}
58+
}
59+
}
60+
for (PCMSolverIndex i = 0; i < nFragments; ++i) {
61+
for (PCMSolverIndex j = 0; j < nSitesPerFragment; ++j) {
62+
Dlambda_(i + nSites, j + i * nSitesPerFragment) = 1.0;
63+
Dlambda_(j + i * nSitesPerFragment, i + nSites) = 1.0;
64+
}
65+
}
66+
TIMER_OFF("Computing Dlambda");
67+
68+
built_ = true;
69+
}
70+
71+
Eigen::VectorXd FQOhno::computeCharge_impl(const Eigen::VectorXd & potential) const {
72+
// We have to add electronegativities on top of the potential
73+
Eigen::VectorXd RHS = Eigen::VectorXd::Zero(Dlambda_.rows());
74+
RHS.head(potential.size()) = potential + mmfq_.chi;
75+
return -Dlambda_.ldlt().solve(RHS).head(potential.size());
76+
}
77+
78+
std::ostream & FQOhno::printSolver(std::ostream & os) {
79+
os << "Fluctuating charge solver type: Ohno";
80+
return os;
81+
}
82+
} // namespace mmfq
83+
} // namespace pcm

‎src/mmfq/FQOhno.hpp‎

Lines changed: 77 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,77 @@
1+
/**
2+
* PCMSolver, an API for the Polarizable Continuum Model
3+
* Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.
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+
#ifndef FQOHNO_HPP
25+
#define FQOHNO_HPP
26+
27+
#include <iosfwd>
28+
29+
#include "Config.hpp"
30+
31+
#include <Eigen/Core>
32+
33+
#include "utils/MMFQ.hpp"
34+
35+
/*! \file FQOhno.hpp
36+
* \class FQOhno
37+
* \brief Solver for MMFQ with Ohno interaction kernel
38+
* \author Roberto Di Remigio
39+
* \date 2017
40+
*
41+
* \note We store the D_\lambda matrix and use a robust Cholesky decomposition
42+
* to solve for the charges.
43+
* This avoids computing and storing the inverse explicitly.
44+
*/
45+
46+
namespace pcm {
47+
namespace mmfq {
48+
class FQOhno __final {
49+
public:
50+
FQOhno() {}
51+
/*! \brief Construct solver
52+
* \param[in] the fluctuating charges force field
53+
*/
54+
FQOhno(const utils::MMFQ & ff);
55+
56+
~FQOhno() {}
57+
friend std::ostream & operator<<(std::ostream & os, FQOhno & solver) {
58+
return solver.printSolver(os);
59+
}
60+
Eigen::VectorXd computeCharge(const Eigen::VectorXd & potential) const {
61+
return computeCharge_impl(potential);
62+
}
63+
64+
private:
65+
bool built_;
66+
utils::MMFQ mmfq_;
67+
/*! D_\lambda matrix, not symmetry blocked */
68+
Eigen::MatrixXd Dlambda_;
69+
70+
void buildSystemMatrix_impl();
71+
Eigen::VectorXd computeCharge_impl(const Eigen::VectorXd & potential) const;
72+
std::ostream & printSolver(std::ostream & os);
73+
};
74+
} // namespace mmfq
75+
} // namespace pcm
76+
77+
#endif // FQOHNO_HPP

‎src/utils/MMFQ.hpp‎

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
/**
2+
* PCMSolver, an API for the Polarizable Continuum Model
3+
* Copyright (C) 2017 Roberto Di Remigio, Luca Frediani and collaborators.
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+
#ifndef MMFQ_HPP
25+
#define MMFQ_HPP
26+
27+
#include <iosfwd>
28+
29+
#include "Config.hpp"
30+
31+
#include <Eigen/Core>
32+
33+
namespace pcm {
34+
namespace utils {
35+
/*! \file MMFQ.hpp
36+
* \struct MMFQ
37+
* \brief POD representing a classical fluctuating charges MM force field
38+
* \author Roberto Di Remigio
39+
* \date 2017
40+
*/
41+
struct MMFQ {
42+
/*! Number of FQ fragments */
43+
PCMSolverIndex nFragments;
44+
/*! Number of FQ sites per MM fragment
45+
* \note This would be three for a standard representation of water
46+
*/
47+
PCMSolverIndex nSitesPerFragment;
48+
/*! FQ electronegativities */
49+
Eigen::VectorXd chi;
50+
/*! FQ hardnesses */
51+
Eigen::VectorXd eta;
52+
/*! FQ sites */
53+
Eigen::Matrix3Xd sites;
54+
};
55+
} // namespace utils
56+
} // namespace pcm
57+
58+
#endif // MMFQ_HPP

‎tests/input/fq.inp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
UNITS = AU
1+
UNITS = Angstrom
22
CAVITY
33
{
44
TYPE = GEPOL

‎tests/mmfq/CMakeLists.txt‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
add_library(mmfq-tests OBJECT mmfq-ohno.cpp)
2+
if(BUILD_CUSTOM_BOOST)
3+
add_dependencies(mmfq-tests custom_boost)
4+
endif()
5+
6+
# Copy reference files to ${PROJECT_BINARY_DIR}/tests/mmfq (aka ${CMAKE_CURRENT_BINARY_DIR})
7+
list(APPEND reference_files FQ-2xH2O.npy)
8+
file(COPY ${reference_files} DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
9+
10+
# mmfq.cpp test
11+
add_Catch_test(mmfq-ohno "mmfq;ohno")
12+

‎tests/mmfq/FQ-2xH2O.npy‎

128 Bytes
Binary file not shown.

0 commit comments

Comments
 (0)