source: ThirdParty/mpqc_open/src/lib/chemistry/molecule/moltest.cc@ 12add2

Action_Thermostats Add_AtomRandomPerturbation Add_RotateAroundBondAction Add_SelectAtomByNameAction Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_StructOpt_integration_tests AutomationFragmentation_failures Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_ChronosMutex Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion GeometryObjects Gui_displays_atomic_force_velocity IndependentFragmentGrids_IntegrationTest JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks RotateToPrincipalAxisSystem_UndoRedo StoppableMakroAction Subpackage_levmar Subpackage_vmg ThirdParty_MPQC_rebuilt_buildsystem TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps Ubuntu_1604_changes stable
Last change on this file since 12add2 was 12add2, checked in by Frederik Heber <heber@…>, 8 years ago

Added missing SRCDIR define to scconfig.h

  • Property mode set to 100644
File size: 10.1 KB
Line 
1//
2// moltest.cc
3//
4// Copyright (C) 1996 Limit Point Systems, Inc.
5//
6// Author: Curtis Janssen <cljanss@limitpt.com>
7// Maintainer: LPS
8//
9// This file is part of the SC Toolkit.
10//
11// The SC Toolkit is free software; you can redistribute it and/or modify
12// it under the terms of the GNU Library General Public License as published by
13// the Free Software Foundation; either version 2, or (at your option)
14// any later version.
15//
16// The SC Toolkit 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 Library General Public License for more details.
20//
21// You should have received a copy of the GNU Library General Public License
22// along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24//
25// The U.S. Government is granted a limited license as per AL 91-7.
26//
27
28#include <scconfig.h>
29
30#include <math.h>
31
32#ifdef HAVE_SSTREAM
33# include <sstream>
34#else
35# include <strstream.h>
36#endif
37
38#include <util/state/stateio.h>
39#include <math/scmat/local.h>
40#include <chemistry/molecule/molecule.h>
41#include <chemistry/molecule/hess.h>
42#include <chemistry/molecule/energy.h>
43#include <chemistry/molecule/coor.h>
44#include <util/state/state_bin.h>
45#include <util/render/object.h>
46#include <util/render/oogl.h>
47#include <util/misc/formio.h>
48#include <chemistry/molecule/formula.h>
49#include <chemistry/molecule/molfreq.h>
50#include <util/group/mstate.h>
51
52using namespace std;
53using namespace sc;
54
55// force linkage
56#include <chemistry/molecule/linkage.h>
57
58void do_displacement(Ref<MolecularCoor>&mc,int i);
59
60void
61print_atominfo(const Ref<AtomInfo> &atominfo,
62 const Ref<AtomInfo> &refatominfo)
63{
64 cout << "Rvdw(H) = " << refatominfo->vdw_radius(1)
65 << " " << atominfo->vdw_radius(1)
66 << "/" << atominfo->vdw_radius_scale()
67 << endl;
68 cout << "Rvdw(C) = " << refatominfo->vdw_radius(6)
69 << " " << atominfo->vdw_radius(6)
70 << "/" << atominfo->vdw_radius_scale()
71 << endl;
72 cout << "Rb(H) = " << refatominfo->bragg_radius(1)
73 << " " << atominfo->bragg_radius(1)
74 << "/" << atominfo->bragg_radius_scale()
75 << endl;
76 cout << "Ra(H) = " << refatominfo->atomic_radius(1)
77 << " " << atominfo->atomic_radius(1)
78 << "/" << atominfo->atomic_radius_scale()
79 << endl;
80 cout << "mass(H) = " << refatominfo->mass(1)
81 << " " << atominfo->mass(1)
82 << endl;
83 cout << "rgb(H) = "
84 << "["
85 << refatominfo->rgb(1,0) << " "
86 << refatominfo->rgb(1,1) << " "
87 << refatominfo->rgb(1,2) << " "
88 << "] ["
89 << atominfo->rgb(1,0) << " "
90 << atominfo->rgb(1,1) << " "
91 << atominfo->rgb(1,2) << " "
92 << "]"
93 << endl;
94}
95
96int
97main(int argc, char **argv)
98{
99 int i;
100
101 // get the message group. first try the commandline and environment
102 Ref<MessageGrp> grp = MessageGrp::initial_messagegrp(argc, argv);
103 if (grp.nonnull())
104 MessageGrp::set_default_messagegrp(grp);
105 else
106 grp = MessageGrp::get_default_messagegrp();
107
108 Ref<KeyVal> kv;
109 if (argc == 2) {
110 kv = new ParsedKeyVal(argv[1]);
111 }
112 else {
113 kv = new ParsedKeyVal(SRCDIR "src/lib/chemistry/molecule/moltest.in");
114 }
115
116 Ref<AtomInfo> atominfo; atominfo << kv->describedclassvalue("atominfo");
117 if (atominfo.nonnull()) {
118 Ref<AtomInfo> refatominfo = new AtomInfo;
119 cout << "-------------- testing atominfo --------------" << endl;
120 if (grp->me() == 0) print_atominfo(atominfo, refatominfo);
121 cout << "saving/restoring atominfo" << endl;
122 StateOutBin so("moltest.1.ckpt");
123 SavableState::save_state(atominfo.pointer(), so);
124 atominfo = 0;
125 so.close();
126 StateInBin si("moltest.1.ckpt");
127 atominfo << SavableState::restore_state(si);
128 if (grp->me() == 0) print_atominfo(atominfo, refatominfo);
129 if (grp->n() > 1) {
130 BcastState b(grp, 0);
131 b.bcast(atominfo);
132 }
133 if (grp->me() == 1) {
134 print_atominfo(atominfo, refatominfo);
135 }
136 }
137
138 Ref<Molecule> mol; mol << kv->describedclassvalue("molecule");
139 if (mol.nonnull()) {
140 cout << "-------------- testing molecule --------------" << endl;
141
142 MolecularFormula formula(mol);
143 cout << "Molecular Formula" << endl << formula.formula() << endl;
144 cout << "Number of Atomtypes" << endl << formula.natomtypes() << endl;
145 cout << "Atomtype, Number of Atoms of This Type" << endl;
146 for(i=0; i<formula.natomtypes(); i++) {
147 cout << formula.Z(i) << "," << formula.nZ(i) << endl;
148 }
149
150 mol->cleanup_molecule();
151 cout << "Clean Molecule:\n";
152 mol->print();
153
154 mol->transform_to_principal_axes();
155 cout << "Clean Molecule wrt principal axes:\n";
156 mol->print();
157
158 int nunique = mol->nunique();
159
160 cout << "nunique=" << nunique << ":";
161 for (i=0; i < nunique; i++) cout << " " << mol->unique(i)+1;
162 cout << endl;
163
164 mol->point_group()->char_table().print();
165
166 cout << "---------- testing molecule save/restore ----------" << endl;
167
168 StateOutBin so("moltest.2.ckpt");
169 cout << "saveing ..." << endl;
170 SavableState::save_state(mol.pointer(),so);
171 mol = 0;
172 so.close();
173 StateInBin si("moltest.2.ckpt");
174 cout << "restoring ..." << endl;
175 mol << SavableState::restore_state(si);
176 cout << "printing restored molecule:" << endl;
177 mol->print();
178 }
179
180 cout << "-------------- initializing render tests --------------" << endl;
181 Ref<Render> ren;
182 ren << kv->describedclassvalue("renderer");
183 Ref<RenderedObject> renmol;
184 renmol << kv->describedclassvalue("renderedmolecule");
185 if (ren.nonnull() && renmol.nonnull()) {
186 cout << "-------------- testing renderer --------------" << endl;
187 ren->render(renmol);
188 }
189
190 //exit(0);
191
192 Ref<SetIntCoor> simp; simp << kv->describedclassvalue("simp");
193 if (simp.nonnull()) {
194 cout << "-------------- testing simp --------------" << endl;
195 Ref<IntCoorGen> gen; gen << kv->describedclassvalue("generator");
196 if (gen.nonnull()) {
197 gen->print();
198 }
199 cout << "simp before update:\n";
200 simp->print_details(mol);
201 simp->update_values(mol);
202 cout << "simp:\n";
203 simp->print_details(mol);
204 }
205
206 // compare the analytic bmatrix to the finite displacement bmatrix
207 Ref<SetIntCoor> bmat_test; bmat_test << kv->describedclassvalue("bmat_test");
208 if (bmat_test.nonnull()) {
209 cout << "-------------- bmat_test --------------" << endl;
210 Ref<SCMatrixKit> kit = SCMatrixKit::default_matrixkit();
211 RefSCDimension dnc(new SCDimension(bmat_test->n()));
212 RefSCDimension dn3(new SCDimension(mol->natom()*3));
213 RefSCMatrix bmatrix(dnc,dn3,kit);
214 RefSCMatrix fd_bmatrix(dnc,dn3,kit);
215 cout << "testing bmat with:\n";
216 bmat_test->update_values(mol);
217 bmat_test->print();
218 bmat_test->bmat(mol,bmatrix);
219 bmat_test->fd_bmat(mol,fd_bmatrix);
220 cout << "test bmatrix:\n";
221 bmatrix.print();
222 cout << "fd bmatrix:\n";
223 fd_bmatrix.print();
224 RefSCMatrix diff = fd_bmatrix - bmatrix;
225 cout << "difference between test and finite displacement bmatrix:\n";
226 diff.print();
227 cout << "% difference between test and finite displacement bmatrix:\n";
228 for (i=0; i<diff.nrow(); i++) {
229 for (int j=0; j<diff.ncol(); j++) {
230 double denom = fabs(fd_bmatrix(i,j));
231 double num = fabs(diff(i,j));
232 if (denom < 0.000001) denom = 0.000001;
233 if (num < 0.00001) diff(i,j) = 0.0;
234 else diff(i,j) = 100.0 * fabs(diff(i,j))/denom;
235 }
236 }
237 diff.print();
238
239 cout << "testing for translational invariance of each coordinate:"
240 << endl;
241 for (i=0; i<bmat_test->n(); i++) {
242 cout << " coor " << scprintf("%2d",i) << ":";
243 for (int j=0; j<3; j++) {
244 double sum = 0.0;
245 for (int k=0; k<mol->natom(); k++) {
246 sum += bmatrix(i,k*3+j);
247 }
248 cout << scprintf(" % 16.12f",sum);
249 }
250 cout << endl;
251 }
252 bmatrix.gi().print("The inverse bmatrix");
253 }
254
255 cout.flush();
256 cerr.flush();
257
258 // now we get ambitious
259 Ref<MolecularCoor> mc; mc << kv->describedclassvalue("molcoor");
260 cout.flush();
261 cerr.flush();
262
263 if (mc.nonnull()) {
264 cout << "-------------- testing molcoor --------------" << endl;
265 mc->print();
266
267 cout.flush();
268 cerr.flush();
269
270 // do_displacement(mc,0);
271 // do_displacement(mc,1);
272 // do_displacement(mc,2);
273 // do_displacement(mc,3);
274
275 Ref<SCMatrixKit> kit = SCMatrixKit::default_matrixkit();
276 RefSymmSCMatrix hessian(mc->dim(),kit);
277 mc->guess_hessian(hessian);
278
279 // cout << "The guess hessian:\n";
280 // hessian.print();
281 }
282
283 Ref<MolecularEnergy> me; me << kv->describedclassvalue("energy");
284 if (me.nonnull()) {
285 cout << "-------------- testing energy --------------" << endl;
286 me->print();
287 }
288
289 Ref<MolecularHessian> molhess; molhess << kv->describedclassvalue("hess");
290 RefSymmSCMatrix xhessian;
291 if (molhess.nonnull()) {
292 xhessian = molhess->cartesian_hessian();
293 }
294
295 Ref<MolecularFrequencies> molfreq;
296 molfreq << kv->describedclassvalue("freq");
297 if (molfreq.nonnull() && xhessian.nonnull()) {
298 cout << "-------------- testing freq --------------" << endl;
299 molfreq->compute_frequencies(xhessian);
300 }
301
302 return 0;
303}
304
305
306void
307do_displacement(Ref<MolecularCoor>&mc,int i)
308{
309 if (i>=mc->dim().n()) return;
310 // now try to displace the geometry
311 RefSCVector internal(mc->dim(),mc->matrixkit());
312 mc->to_internal(internal);
313 cout << "The initial internal coordinates:\n";
314 internal.print();
315 internal(i) = internal(i) + 0.2;
316 cout << "The new internal coordinates:\n";
317 internal.print();
318 mc->to_cartesian(internal);
319 mc->to_internal(internal);
320 cout << "The actual new internal coordinates:\n";
321 internal.print();
322}
323
324/////////////////////////////////////////////////////////////////////////////
325
326// Local Variables:
327// mode: c++
328// c-file-style: "CLJ"
329// End:
Note: See TracBrowser for help on using the repository browser.