ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/Thermo.cpp
Revision: 218
Committed: Sun Dec 29 19:11:05 2002 UTC (21 years, 6 months ago) by chuckv
File size: 5728 byte(s)
Log Message:
Changed thermo getKinetic and getPotential for MPI.
Added header for SPRNG random number generator class.

File Contents

# Content
1 #include <cmath>
2 #include <mpi++.h>
3
4 #include "Thermo.hpp"
5 #include "SRI.hpp"
6 #include "LRI.hpp"
7 #include "Integrator.hpp"
8
9
10 double Thermo::getKinetic(){
11
12 const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
13 double vx2, vy2, vz2;
14 double kinetic, v_sqr;
15 int kl;
16 double jx2, jy2, jz2; // the square of the angular momentums
17
18 DirectionalAtom *dAtom;
19
20 int n_atoms;
21 double kinetic_global;
22 Atom** atoms;
23
24
25 n_atoms = entry_plug->n_atoms;
26 atoms = entry_plug->atoms;
27
28 kinetic = 0.0;
29 kinetic_global = 0.0;
30 for( kl=0; kl < n_atoms; kl++ ){
31
32 vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
33 vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
34 vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();
35
36 v_sqr = vx2 + vy2 + vz2;
37 kinetic += atoms[kl]->getMass() * v_sqr;
38
39 if( atoms[kl]->isDirectional() ){
40
41 dAtom = (DirectionalAtom *)atoms[kl];
42
43 jx2 = dAtom->getJx() * dAtom->getJx();
44 jy2 = dAtom->getJy() * dAtom->getJy();
45 jz2 = dAtom->getJz() * dAtom->getJz();
46
47 kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
48 + (jz2 / dAtom->getIzz());
49 }
50 }
51 #ifdef IS_MPI
52 MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
53 kinetic = kinetic_global;
54 #endif
55
56 kinetic = kinetic * 0.5 / e_convert;
57
58 return kinetic;
59 }
60
61 double Thermo::getPotential(){
62
63 double potential;
64 double potential_global;
65 int el, nSRI;
66 SRI** sris;
67
68 sris = entry_plug->sr_interactions;
69 nSRI = entry_plug->n_SRI;
70
71 potential = 0.0;
72 potential_global = 0.0;
73 potential += entry_plug->longRange->get_potential();;
74
75 // std::cerr << "long range potential: " << potential << "\n";
76 for( el=0; el<nSRI; el++ ){
77
78 potential += sris[el]->get_potential();
79 }
80
81 // Get total potential for entire system from MPI.
82 #ifdef IS_MPI
83 MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
84 potential = potential_global;
85 #endif
86
87 return potential;
88 }
89
90 double Thermo::getTotalE(){
91
92 double total;
93
94 total = this->getKinetic() + this->getPotential();
95 return total;
96 }
97
98 double Thermo::getTemperature(){
99
100 const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
101 double temperature;
102
103 int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
104 - entry_plug->n_constraints - 3;
105
106 temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
107 return temperature;
108 }
109
110 double Thermo::getPressure(){
111
112 // const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
113 // const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
114 // const double conv_A_m = 1.0E-10; //convert A -> m
115
116 return 0.0;
117 }
118
119 void Thermo::velocitize() {
120
121 double x,y;
122 double vx, vy, vz;
123 double jx, jy, jz;
124 int i, vr, vd; // velocity randomizer loop counters
125 double vdrift[3];
126 double mtot = 0.0;
127 double vbar;
128 const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
129 double av2;
130 double kebar;
131 int ndf; // number of degrees of freedom
132 int ndfRaw; // the raw number of degrees of freedom
133 int n_atoms;
134 Atom** atoms;
135 DirectionalAtom* dAtom;
136 double temperature;
137 int n_oriented;
138 int n_constraints;
139
140 atoms = entry_plug->atoms;
141 n_atoms = entry_plug->n_atoms;
142 temperature = entry_plug->target_temp;
143 n_oriented = entry_plug->n_oriented;
144 n_constraints = entry_plug->n_constraints;
145
146
147 ndfRaw = 3 * n_atoms + 3 * n_oriented;
148 ndf = ndfRaw - n_constraints - 3;
149 kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
150
151 for(vr = 0; vr < n_atoms; vr++){
152
153 // uses equipartition theory to solve for vbar in angstrom/fs
154
155 av2 = 2.0 * kebar / atoms[vr]->getMass();
156 vbar = sqrt( av2 );
157
158 // vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
159
160 // picks random velocities from a gaussian distribution
161 // centered on vbar
162
163 x = drand48();
164 y = drand48();
165 vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
166
167 x = drand48();
168 y = drand48();
169 vy = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
170
171 x = drand48();
172 y = drand48();
173 vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
174
175 atoms[vr]->set_vx( vx );
176 atoms[vr]->set_vy( vy );
177 atoms[vr]->set_vz( vz );
178 }
179
180 // Corrects for the center of mass drift.
181 // sums all the momentum and divides by total mass.
182
183 mtot = 0.0;
184 vdrift[0] = 0.0;
185 vdrift[1] = 0.0;
186 vdrift[2] = 0.0;
187 for(vd = 0; vd < n_atoms; vd++){
188
189 vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
190 vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
191 vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
192
193 mtot = mtot + atoms[vd]->getMass();
194 }
195
196 for (vd = 0; vd < 3; vd++) {
197 vdrift[vd] = vdrift[vd] / mtot;
198 }
199
200 for(vd = 0; vd < n_atoms; vd++){
201
202 vx = atoms[vd]->get_vx();
203 vy = atoms[vd]->get_vy();
204 vz = atoms[vd]->get_vz();
205
206
207 vx -= vdrift[0];
208 vy -= vdrift[1];
209 vz -= vdrift[2];
210
211 atoms[vd]->set_vx(vx);
212 atoms[vd]->set_vy(vy);
213 atoms[vd]->set_vz(vz);
214 }
215 if( n_oriented ){
216
217 for( i=0; i<n_atoms; i++ ){
218
219 if( atoms[i]->isDirectional() ){
220
221 dAtom = (DirectionalAtom *)atoms[i];
222
223 vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
224 x = drand48();
225 y = drand48();
226 jx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
227
228 vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
229 x = drand48();
230 y = drand48();
231 jy = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
232
233 vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
234 x = drand48();
235 y = drand48();
236 jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
237
238 dAtom->setJx( jx );
239 dAtom->setJy( jy );
240 dAtom->setJz( jz );
241 }
242 }
243 }
244 }