1 /**********************************************************************
2 data_utilities.cpp - Global data and resource file parsers.
3
4 Copyright (C) 1998-2001 by OpenEye Scientific Software, Inc.
5 Copyright (C) 2015 by David van der Spoel
6
7 This file is part of the Open Babel project.
8 For more information, see <http://openbabel.org/>
9
10 This program is free software; you can redistribute it and/or modify
11 it under the terms of the GNU General Public License as published by
12 the Free Software Foundation version 2 of the License.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18 ***********************************************************************/
19
20 #ifdef WIN32
21 #pragma warning (disable : 4786)
22 #endif
23 #include <cstdlib>
24 #include <cstring>
25 #include <fstream>
26 #include <iostream>
27 #include <openbabel/babelconfig.h>
28 #include <openbabel/data_utilities.h>
29 #include <openbabel/mol.h>
30 #include <openbabel/generic.h>
31 #include <openbabel/locale.h>
32
33 namespace OpenBabel {
34
extract_thermochemistry(OpenBabel::OBMol & mol,bool bVerbose,int * Nsymm,int Nrotbonds,double dBdT,double * temperature,double * DeltaHf0,double * DeltaHfT,double * DeltaGfT,double * DeltaSfT,double * S0T,double * CVT,double * CPT,std::vector<double> & Scomponents,double * ZPVE)35 bool extract_thermochemistry(OpenBabel::OBMol &mol,
36 bool bVerbose,
37 int *Nsymm,
38 int Nrotbonds,
39 double dBdT,
40 double *temperature,
41 double *DeltaHf0,
42 double *DeltaHfT,
43 double *DeltaGfT,
44 double *DeltaSfT,
45 double *S0T,
46 double *CVT,
47 double *CPT,
48 std::vector<double> &Scomponents,
49 double *ZPVE)
50 {
51 enum kkTYPE {kkDH, kkDG, kkDS, kkS0, kkCV, kkSt, kkSr, kkSv, kkZP};
52 typedef struct {
53 std::string term;
54 kkTYPE kk;
55 } energy_unit;
56 double St = 0, Sr = 0, Sv = 0, Sconf = 0, Ssymm = 0;
57 double Rgas = 1.9872041;
58 int RotSymNum = 1;
59 OpenBabel::OBRotationData* rd;
60
61 rd = (OpenBabel::OBRotationData*)mol.GetData("RotationData");
62 if (nullptr != rd)
63 {
64 RotSymNum = rd->GetSymmetryNumber();
65 if (bVerbose)
66 {
67 printf("Found symmetry number %d in input file.\n", RotSymNum);
68 }
69 }
70 else if (bVerbose)
71 {
72 printf("Using default symmetry number %d\n", RotSymNum);
73 }
74 if ((*Nsymm > 0) && (*Nsymm != RotSymNum))
75 {
76 // Rgas in cal/mol K http://en.wikipedia.org/wiki/Gas_constant
77 Ssymm = -Rgas*log((1.0* *Nsymm)/RotSymNum);
78 RotSymNum = *Nsymm;
79 if (bVerbose)
80 {
81 printf("Changing symmetry number to %d\n", RotSymNum);
82 }
83 }
84 else if (*Nsymm == 0)
85 {
86 *Nsymm = RotSymNum;
87 }
88 if (Nrotbonds > 0)
89 {
90 Sconf = Rgas*Nrotbonds*log(3.0);
91 }
92 energy_unit eu[] = {
93 { "zpe", kkZP },
94 { "DeltaHform", kkDH },
95 { "DeltaGform", kkDG },
96 { "DeltaSform", kkDS },
97 { "S0", kkS0 },
98 { "cv", kkCV },
99 { "Strans", kkSt },
100 { "Srot", kkSr },
101 { "Svib", kkSv }
102 };
103 #define NEU (sizeof(eu)/sizeof(eu[0]))
104 int found = 0;
105 std::vector<OpenBabel::OBGenericData*> obdata = mol.GetData();
106 for(std::vector<OpenBabel::OBGenericData*>::iterator j = obdata.begin(); (j<obdata.end()); ++j)
107 {
108 std::string term = (*j)->GetAttribute();
109 double value = atof((*j)->GetValue().c_str());
110 double T = 0;
111 {
112 size_t lh = term.find("(");
113 size_t rh = term.find("K)");
114 double TT = atof(term.substr(lh+1,rh-lh-1).c_str());
115 if (0 != TT)
116 {
117 if (0 == T)
118 {
119 T = TT;
120 *temperature = TT;
121 }
122 else
123 {
124 std::cerr << "Different T in the input file, found " << T << " before and now " << TT << ". Output maybe inconsistent." << std::endl;
125 T = TT;
126 }
127 }
128 }
129 for(unsigned int i = 0; (i<NEU); i++)
130 {
131 if (strstr(term.c_str(), eu[i].term.c_str()) != nullptr)
132 {
133 switch (eu[i].kk)
134 {
135 case kkZP:
136 {
137 *ZPVE = value;
138 }
139 break;
140 case kkDH:
141 if (0 == T)
142 {
143 *DeltaHf0 = value;
144 }
145 else
146 {
147 *DeltaHfT = value;
148 }
149 found ++;
150 break;
151 case kkDG:
152 *DeltaGfT = value - T*(Ssymm+Sconf)/1000;
153 found ++;
154 break;
155 case kkDS:
156 *DeltaSfT = value + Ssymm + Sconf;
157 found ++;
158 break;
159 case kkS0:
160 *S0T = value + Ssymm + Sconf;
161 found ++;
162 break;
163 case kkSt:
164 St = value;
165 found ++;
166 break;
167 case kkSr:
168 Sr = value;
169 found ++;
170 break;
171 case kkSv:
172 Sv = value;
173 found ++;
174 break;
175 case kkCV:
176 *CVT = value;
177 found++;
178 break;
179 default:
180 break;
181 }
182 }
183 }
184 }
185 double P = 16.605/4.184; // Convert pressure to kcal/mol
186 *CPT = *CVT + Rgas + (2*P*dBdT + pow(P*dBdT, 2.0)/Rgas);
187
188 Scomponents.push_back(St);
189 Scomponents.push_back(Sr);
190 Scomponents.push_back(Sv);
191 Scomponents.push_back(Ssymm);
192 Scomponents.push_back(Sconf);
193 if (bVerbose && (Ssymm != 0))
194 {
195 printf("Applyied symmetry correction to free energy of %g kcal/mol\n",
196 -(*temperature*Ssymm)/1000);
197 }
198 if (bVerbose && (Sconf != 0))
199 {
200 printf("Applyied conformational correction to free energy of %g kcal/mol\n",
201 -(*temperature*Sconf)/1000);
202 }
203 return (found == 9);
204 }
205
206 }
207
208 //! \file data_utilities.cpp
209 //! \brief Data related tools and utilities
210