1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: qmmm.cc
4 // Copyright (C) 2017 Toru Shiozaki
5 //
6 // Author: Jae Woo Park <jwpk1201@northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program. If not, see <http://www.gnu.org/licenses/>.
23 //
24
25
26 #include <functional>
27 #include <typeinfo>
28 #include <fstream>
29 #include <string>
30 #include <algorithm>
31 #include <src/grad/gradeval.h>
32 #include <src/util/timer.h>
33 #include <src/opt/opt.h>
34
35 using namespace std;
36 using namespace bagel;
37
edit_tinker_input(shared_ptr<const Geometry> current) const38 void QMMM_Tinker::edit_tinker_input(shared_ptr<const Geometry> current) const {
39 system("mv -f tinkin.xyz tinkin.xyz.old");
40 ifstream fs_tinker_qmmm_old("tinkin.xyz.old");
41 ofstream fs_tinker_qmmm("tinkin.xyz");
42 string line;
43 getline(fs_tinker_qmmm_old, line);
44 fs_tinker_qmmm << line << endl;
45 {
46 int i = 0;
47 while (getline(fs_tinker_qmmm_old, line)) {
48 stringstream ss(line);
49 int dum;
50 string atomnm;
51 double x,y,z;
52 ss >> dum >> atomnm >> x >> y >> z;
53 fs_tinker_qmmm << setw(6) << dum << setw(3) << atomnm << setw(20) << setprecision(10) << current->xyz()->element(0, i) * au2angstrom__ <<
54 setw(20) << setprecision(10) << current->xyz()->element(1, i) * au2angstrom__ <<
55 setw(20) << setprecision(10) << current->xyz()->element(2, i) * au2angstrom__;
56 // get the rest dummies and print it to TINKER input
57 while (ss >> dum)
58 fs_tinker_qmmm << setw(6) << dum;
59 fs_tinker_qmmm << endl;
60 ++i;
61 }
62 }
63 }
64
65
edit_input(shared_ptr<const Geometry> current) const66 void QMMM_Tinker::edit_input(shared_ptr<const Geometry> current) const {
67 // generate TINKER input & run TINKER testgrad.x
68 // (tinker inputs are prepared in tinker1 and tinker2 directories)
69
70 Timer timer;
71 int natom = current->natom();
72
73 chdir("tinker1");
74
75 {
76 edit_tinker_input(current);
77 system("testgrad -k tinkin.key tinkin.xyz y n n > gradientls");
78 system("grep 'Total Potential' gradientls | cut -c 29- > ../energy_qmmm");
79 stringstream ss;
80 ss << "grep -A " << natom + 3 << " 'Cartesian Gradient Breakdown over Individual Atoms' gradientls | grep -v 'Gradient' | grep 'Anlyt' | cut -c 8- > ../grad_qmmm";
81 string pp = ss.str();
82 system(pp.c_str());
83 }
84
85 chdir("../");
86 timer.tick_print("Running TINKER for whole region");
87 chdir("tinker2");
88
89 {
90 edit_tinker_input(current);
91 system("testgrad -k tinkin.key tinkin.xyz y n n > gradientls");
92 system("grep 'Total Potential' gradientls | cut -c 29- > ../energy_qm");
93 stringstream ss;
94 ss << "grep -A " << natom + 3 << " 'Cartesian Gradient Breakdown over Individual Atoms' gradientls | grep -v 'Gradient' | grep 'Anlyt' | cut -c 8- > ../grad_qm";
95 string pp = ss.str();
96 system(pp.c_str());
97 }
98
99 chdir("../");
100 timer.tick_print("Running TINKER for QM region");
101 }
102
103
do_grad(const int natom) const104 tuple<double,shared_ptr<GradFile>> QMMM_Tinker::do_grad(const int natom) const {
105 Timer timer;
106 double mmen;
107 auto out = make_shared<GradFile>(natom);
108
109 // parse all the outputs
110 double mmen_1;
111 {
112 ifstream fs_energy_qmmm("energy_qmmm");
113 if (!fs_energy_qmmm.is_open())
114 throw runtime_error("energy_qmmm cannot be opened");
115 string line;
116 {
117 getline(fs_energy_qmmm, line);
118 stringstream ss(line);
119 ss >> mmen_1;
120 }
121 }
122
123 double mmen_2;
124 {
125 ifstream fs_energy_qm("energy_qm");
126 if (!fs_energy_qm.is_open())
127 throw runtime_error("energy_qm cannot be opened");
128 string line;
129 {
130 getline(fs_energy_qm, line);
131 stringstream ss(line);
132 ss >> mmen_2;
133 }
134 }
135
136 auto grad_1 = make_shared<GradFile>(natom);
137 {
138 ifstream fs_grad_qmmm("grad_qmmm");
139 if (!fs_grad_qmmm.is_open())
140 throw runtime_error("grad_qmmm cannot be opened");
141 string line;
142 {
143 int i = 0;
144 while (getline(fs_grad_qmmm, line)) {
145 int dum;
146 stringstream ss(line);
147 ss >> dum >> grad_1->element(0, i) >> grad_1->element(1, i) >> grad_1->element(2, i);
148 ++i;
149 }
150 }
151 }
152
153 auto grad_2 = make_shared<GradFile>(natom);
154 {
155 ifstream fs_grad_qm("grad_qm");
156 if (!fs_grad_qm.is_open())
157 throw runtime_error("grad_qm cannot be opened");
158 string line;
159 {
160 int i = 0;
161 while (getline(fs_grad_qm, line)) {
162 int dum;
163 stringstream ss(line);
164 ss >> dum >> grad_2->element(0, i) >> grad_2->element(1, i) >> grad_2->element(2, i);
165 ++i;
166 }
167 }
168 }
169
170 // energy : kcal / mol
171 mmen = (mmen_1 - mmen_2) * kcal2kj__ / (au2kjmol__);
172
173 stringstream ss; ss << "MM energy = " << setw(10) << setprecision(5) << mmen;
174 timer.tick_print(ss.str());
175
176 // gradient : kcal / mol / angstrom
177 grad_1->scale(kcal2kj__ * au2angstrom__ / au2kjmol__);
178 grad_2->scale(kcal2kj__ * au2angstrom__ / au2kjmol__);
179
180 *out = *grad_1 - *grad_2;
181
182 // return the energy and gradient
183 return tie(mmen, out);
184 }
185