ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/NPT.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/NPT.cpp (file contents):
Revision 857 by mmeineke, Fri Nov 7 17:09:48 2003 UTC vs.
Revision 1268 by tim, Fri Jun 11 17:16:21 2004 UTC

# Line 62 | Line 62 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
62      integralOfChidt = integralOfChidtValue->getData();
63    }
64  
65 <  oldPos = new double[3*nAtoms];
66 <  oldVel = new double[3*nAtoms];
67 <  oldJi = new double[3*nAtoms];
68 < #ifdef IS_MPI
69 <  Nparticles = mpiSim->getTotAtoms();
70 < #else
71 <  Nparticles = theInfo->n_atoms;
72 < #endif
65 >  oldPos = new double[3*integrableObjects.size()];
66 >  oldVel = new double[3*integrableObjects.size()];
67 >  oldJi = new double[3*integrableObjects.size()];
68  
69   }
70  
# Line 83 | Line 78 | template<typename T> void NPT<T>::moveA() {
78  
79    //new version of NPT
80    int i, j, k;
86 –  DirectionalAtom* dAtom;
81    double Tb[3], ji[3];
82    double mass;
83    double vel[3], pos[3], frc[3];
# Line 101 | Line 95 | template<typename T> void NPT<T>::moveA() {
95  
96    calcVelScale();
97    
98 <  for( i=0; i<nAtoms; i++ ){
98 >  for( i=0; i<integrableObjects.size(); i++ ){
99  
100 <    atoms[i]->getVel( vel );
101 <    atoms[i]->getFrc( frc );
100 >    integrableObjects[i]->getVel( vel );
101 >    integrableObjects[i]->getFrc( frc );
102  
103 <    mass = atoms[i]->getMass();
103 >    mass = integrableObjects[i]->getMass();
104  
105      getVelScaleA( sc, vel );
106  
# Line 117 | Line 111 | template<typename T> void NPT<T>::moveA() {
111  
112      }
113  
114 <    atoms[i]->setVel( vel );
114 >    integrableObjects[i]->setVel( vel );
115  
116 <    if( atoms[i]->isDirectional() ){
116 >    if( integrableObjects[i]->isDirectional() ){
117  
124 –      dAtom = (DirectionalAtom *)atoms[i];
125 –
118        // get and convert the torque to body frame
119  
120 <      dAtom->getTrq( Tb );
121 <      dAtom->lab2Body( Tb );
120 >      integrableObjects[i]->getTrq( Tb );
121 >      integrableObjects[i]->lab2Body( Tb );
122  
123        // get the angular momentum, and propagate a half step
124  
125 <      dAtom->getJ( ji );
125 >      integrableObjects[i]->getJ( ji );
126  
127        for (j=0; j < 3; j++)
128          ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
129  
130 <      this->rotationPropagation( dAtom, ji );
130 >      this->rotationPropagation( integrableObjects[i], ji );
131  
132 <      dAtom->setJ( ji );
132 >      integrableObjects[i]->setJ( ji );
133      }
134    }
135  
# Line 150 | Line 142 | template<typename T> void NPT<T>::moveA() {
142    integralOfChidt += dt2*chi;
143  
144    //save the old positions
145 <  for(i = 0; i < nAtoms; i++){
146 <    atoms[i]->getPos(pos);
145 >  for(i = 0; i < integrableObjects.size(); i++){
146 >    integrableObjects[i]->getPos(pos);
147      for(j = 0; j < 3; j++)
148        oldPos[i*3 + j] = pos[j];
149    }
# Line 160 | Line 152 | template<typename T> void NPT<T>::moveA() {
152  
153    for(k = 0; k < 5; k ++){
154  
155 <    for(i =0 ; i < nAtoms; i++){
155 >    for(i =0 ; i < integrableObjects.size(); i++){
156  
157 <      atoms[i]->getVel(vel);
158 <      atoms[i]->getPos(pos);
157 >      integrableObjects[i]->getVel(vel);
158 >      integrableObjects[i]->getPos(pos);
159  
160        this->getPosScale( pos, COM, i, sc );
161  
162        for(j = 0; j < 3; j++)
163          pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
164  
165 <      atoms[i]->setPos( pos );
165 >      integrableObjects[i]->setPos( pos );
166      }
167  
168 <    if (nConstrained){
177 <      constrainA();
178 <    }
168 >    consFramework->doConstrainA();
169    }
170  
171  
# Line 188 | Line 178 | template<typename T> void NPT<T>::moveB( void ){
178  
179    //new version of NPT
180    int i, j, k;
191 –  DirectionalAtom* dAtom;
181    double Tb[3], ji[3], sc[3];
182    double vel[3], frc[3];
183    double mass;
184  
185    // Set things up for the iteration:
186  
187 <  for( i=0; i<nAtoms; i++ ){
187 >  for( i=0; i<integrableObjects.size(); i++ ){
188  
189 <    atoms[i]->getVel( vel );
189 >    integrableObjects[i]->getVel( vel );
190  
191      for (j=0; j < 3; j++)
192        oldVel[3*i + j]  = vel[j];
193  
194 <    if( atoms[i]->isDirectional() ){
194 >    if( integrableObjects[i]->isDirectional() ){
195  
196 <      dAtom = (DirectionalAtom *)atoms[i];
196 >      integrableObjects[i]->getJ( ji );
197  
209 –      dAtom->getJ( ji );
210 –
198        for (j=0; j < 3; j++)
199          oldJi[3*i + j] = ji[j];
200  
# Line 229 | Line 216 | template<typename T> void NPT<T>::moveB( void ){
216      this->evolveEtaB();
217      this->calcVelScale();
218  
219 <    for( i=0; i<nAtoms; i++ ){
219 >    for( i=0; i<integrableObjects.size(); i++ ){
220  
221 <      atoms[i]->getFrc( frc );
222 <      atoms[i]->getVel(vel);
221 >      integrableObjects[i]->getFrc( frc );
222 >      integrableObjects[i]->getVel(vel);
223  
224 <      mass = atoms[i]->getMass();
224 >      mass = integrableObjects[i]->getMass();
225  
226        getVelScaleB( sc, i );
227  
# Line 242 | Line 229 | template<typename T> void NPT<T>::moveB( void ){
229        for (j=0; j < 3; j++)
230          vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
231  
232 <      atoms[i]->setVel( vel );
232 >      integrableObjects[i]->setVel( vel );
233  
234 <      if( atoms[i]->isDirectional() ){
234 >      if( integrableObjects[i]->isDirectional() ){
235  
249 –        dAtom = (DirectionalAtom *)atoms[i];
250 –
236          // get and convert the torque to body frame
237  
238 <        dAtom->getTrq( Tb );
239 <        dAtom->lab2Body( Tb );
238 >        integrableObjects[i]->getTrq( Tb );
239 >        integrableObjects[i]->lab2Body( Tb );
240  
241          for (j=0; j < 3; j++)
242            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
243  
244 <          dAtom->setJ( ji );
244 >          integrableObjects[i]->setJ( ji );
245        }
246      }
247  
248 <    if (nConstrained){
264 <      constrainB();
265 <    }
248 >    consFramework->doConstrainA();
249  
250      if ( this->chiConverged() && this->etaConverged() ) break;
251    }
# Line 338 | Line 321 | template<typename T> int NPT<T>::readyCheck() {
321  
322    if (!have_tau_barostat) {
323      sprintf( painCave.errMsg,
324 <             "NPT error: If you use the NPT\n"
325 <             "   integrator, you must set tauBarostat.\n");
324 >             "If you use the NPT integrator, you must set tauBarostat.\n");
325 >    painCave.severity = OOPSE_ERROR;
326      painCave.isFatal = 1;
327      simError();
328      return -1;
# Line 347 | Line 330 | template<typename T> int NPT<T>::readyCheck() {
330  
331    if (!have_chi_tolerance) {
332      sprintf( painCave.errMsg,
333 <             "NPT warning: setting chi tolerance to 1e-6\n");
333 >             "Setting chi tolerance to 1e-6 in NPT integrator\n");
334      chiTolerance = 1e-6;
335      have_chi_tolerance = 1;
336 +    painCave.severity = OOPSE_INFO;
337      painCave.isFatal = 0;
338      simError();
339    }
340  
341    if (!have_eta_tolerance) {
342      sprintf( painCave.errMsg,
343 <             "NPT warning: setting eta tolerance to 1e-6\n");
343 >             "Setting eta tolerance to 1e-6 in NPT integrator");
344      etaTolerance = 1e-6;
345      have_eta_tolerance = 1;
346 +    painCave.severity = OOPSE_INFO;
347      painCave.isFatal = 0;
348      simError();
349    }
350  
351    // We need NkBT a lot, so just set it here: This is the RAW number
352 <  // of particles, so no subtraction or addition of constraints or
352 >  // of integrableObjects, so no subtraction or addition of constraints or
353    // orientational degrees of freedom:
354  
355 <  NkBT = (double)Nparticles * kB * targetTemp;
355 >  NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp;
356  
357    // fkBT is used because the thermostat operates on more degrees of freedom
358    // than the barostat (when there are particles with orientational degrees
359 <  // of freedom).  ndf = 3 * (n_atoms + n_oriented -1) - n_constraint - nZcons
359 >  // of freedom).  
360  
361 <  fkBT = (double)info->ndf * kB * targetTemp;
361 >  fkBT = (double)(info->getNDF()) * kB * targetTemp;
362  
363    tt2 = tauThermostat * tauThermostat;
364    tb2 = tauBarostat * tauBarostat;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines