# | Line 224 | Line 224 | void Integrator::integrate( void ){ | |
---|---|---|
224 | pos = Atom::getPosArray(); | |
225 | vel = Atom::getVelArray(); | |
226 | frc = Atom::getFrcArray(); | |
227 | – | trq = Atom::getTrqArray(); |
228 | – | Amat = Atom::getAmatArray(); |
227 | ||
228 | while( currTime < runTime ){ | |
229 | ||
# | Line 234 | Line 232 | void Integrator::integrate( void ){ | |
232 | calcStress = 1; | |
233 | } | |
234 | ||
235 | < | std::cerr << "calcPot = " << calcPot << "; calcStress = " |
238 | < | << calcStress << "\n"; |
235 | > | std::cerr << currTime << "\n"; |
236 | ||
237 | integrateStep( calcPot, calcStress ); | |
238 | ||
# | Line 282 | Line 279 | void Integrator::integrateStep( int calcPot, int calcS | |
279 | ||
280 | preMove(); | |
281 | moveA(); | |
282 | < | if( nConstrained ) constrainA(); |
282 | > | //if( nConstrained ) constrainA(); |
283 | ||
284 | // calc forces | |
285 | ||
# | Line 304 | Line 301 | void Integrator::moveA( void ){ | |
301 | double Tb[3]; | |
302 | double ji[3]; | |
303 | double angle; | |
304 | < | double A[3][3]; |
304 | > | double A[3][3], At[3][3]; |
305 | ||
306 | ||
307 | for( i=0; i<nAtoms; i++ ){ | |
# | Line 315 | Line 312 | void Integrator::moveA( void ){ | |
312 | for( j=atomIndex; j<(atomIndex+3); j++ ) | |
313 | vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert; | |
314 | ||
318 | – | std::cerr<< "MoveA vel[" << i << "] = " |
319 | – | << vel[atomIndex] << "\t" |
320 | – | << vel[atomIndex+1]<< "\t" |
321 | – | << vel[atomIndex+2]<< "\n"; |
315 | ||
316 | // position whole step | |
317 | for( j=atomIndex; j<(atomIndex+3); j++ ) pos[j] += dt * vel[j]; | |
318 | ||
319 | ||
327 | – | std::cerr<< "MoveA pos[" << i << "] = " |
328 | – | << pos[atomIndex] << "\t" |
329 | – | << pos[atomIndex+1]<< "\t" |
330 | – | << pos[atomIndex+2]<< "\n"; |
331 | – | |
320 | if( atoms[i]->isDirectional() ){ | |
321 | ||
322 | dAtom = (DirectionalAtom *)atoms[i]; | |
# | Line 338 | Line 326 | void Integrator::moveA( void ){ | |
326 | Tb[0] = dAtom->getTx(); | |
327 | Tb[1] = dAtom->getTy(); | |
328 | Tb[2] = dAtom->getTz(); | |
329 | < | |
329 | > | |
330 | dAtom->lab2Body( Tb ); | |
331 | < | |
331 | > | |
332 | // get the angular momentum, and propagate a half step | |
333 | ||
334 | ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * eConvert; | |
# | Line 349 | Line 337 | void Integrator::moveA( void ){ | |
337 | ||
338 | // use the angular velocities to propagate the rotation matrix a | |
339 | // full time step | |
352 | – | |
353 | – | // get the atom's rotation matrix |
354 | – | |
355 | – | A[0][0] = dAtom->getAxx(); |
356 | – | A[0][1] = dAtom->getAxy(); |
357 | – | A[0][2] = dAtom->getAxz(); |
340 | ||
359 | – | A[1][0] = dAtom->getAyx(); |
360 | – | A[1][1] = dAtom->getAyy(); |
361 | – | A[1][2] = dAtom->getAyz(); |
362 | – | |
363 | – | A[2][0] = dAtom->getAzx(); |
364 | – | A[2][1] = dAtom->getAzy(); |
365 | – | A[2][2] = dAtom->getAzz(); |
366 | – | |
341 | // rotate about the x-axis | |
342 | angle = dt2 * ji[0] / dAtom->getIxx(); | |
343 | < | this->rotate( 1, 2, angle, ji, A ); |
344 | < | |
343 | > | this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); |
344 | > | |
345 | // rotate about the y-axis | |
346 | angle = dt2 * ji[1] / dAtom->getIyy(); | |
347 | < | this->rotate( 2, 0, angle, ji, A ); |
347 | > | this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); |
348 | ||
349 | // rotate about the z-axis | |
350 | angle = dt * ji[2] / dAtom->getIzz(); | |
351 | < | this->rotate( 0, 1, angle, ji, A ); |
351 | > | this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] ); |
352 | ||
353 | // rotate about the y-axis | |
354 | angle = dt2 * ji[1] / dAtom->getIyy(); | |
355 | < | this->rotate( 2, 0, angle, ji, A ); |
355 | > | this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] ); |
356 | ||
357 | // rotate about the x-axis | |
358 | angle = dt2 * ji[0] / dAtom->getIxx(); | |
359 | < | this->rotate( 1, 2, angle, ji, A ); |
359 | > | this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); |
360 | ||
361 | dAtom->setJx( ji[0] ); | |
362 | dAtom->setJy( ji[1] ); | |
363 | dAtom->setJz( ji[2] ); | |
364 | + | |
365 | + | std::cerr << "Amat[" << i << "]\n"; |
366 | + | info->printMat9( &Amat[aMatIndex] ); |
367 | + | |
368 | + | std::cerr << "ji[" << i << "]\t" |
369 | + | << ji[0] << "\t" |
370 | + | << ji[1] << "\t" |
371 | + | << ji[2] << "\n"; |
372 | + | |
373 | } | |
374 | ||
375 | } | |
# | Line 395 | Line 378 | void Integrator::moveB( void ){ | |
378 | ||
379 | void Integrator::moveB( void ){ | |
380 | int i,j,k; | |
381 | < | int atomIndex; |
381 | > | int atomIndex, aMatIndex; |
382 | DirectionalAtom* dAtom; | |
383 | double Tb[3]; | |
384 | double ji[3]; | |
385 | ||
386 | for( i=0; i<nAtoms; i++ ){ | |
387 | atomIndex = i * 3; | |
388 | + | aMatIndex = i * 9; |
389 | ||
390 | // velocity half step | |
391 | for( j=atomIndex; j<(atomIndex+3); j++ ) | |
392 | vel[j] += ( dt2 * frc[j] / atoms[i]->getMass() ) * eConvert; | |
393 | ||
394 | < | std::cerr<< "MoveB vel[" << i << "] = " |
411 | < | << vel[atomIndex] << "\t" |
412 | < | << vel[atomIndex+1]<< "\t" |
413 | < | << vel[atomIndex+2]<< "\n"; |
414 | < | |
415 | < | |
394 | > | |
395 | if( atoms[i]->isDirectional() ){ | |
396 | ||
397 | dAtom = (DirectionalAtom *)atoms[i]; | |
# | Line 423 | Line 402 | void Integrator::moveB( void ){ | |
402 | Tb[1] = dAtom->getTy(); | |
403 | Tb[2] = dAtom->getTz(); | |
404 | ||
405 | + | std::cerr << "TrqB[" << i << "]\t" |
406 | + | << Tb[0] << "\t" |
407 | + | << Tb[1] << "\t" |
408 | + | << Tb[2] << "\n"; |
409 | + | |
410 | dAtom->lab2Body( Tb ); | |
411 | ||
412 | // get the angular momentum, and complete the angular momentum | |
# | Line 435 | Line 419 | void Integrator::moveB( void ){ | |
419 | dAtom->setJx( ji[0] ); | |
420 | dAtom->setJy( ji[1] ); | |
421 | dAtom->setJz( ji[2] ); | |
422 | + | |
423 | + | |
424 | + | std::cerr << "Amat[" << i << "]\n"; |
425 | + | info->printMat9( &Amat[aMatIndex] ); |
426 | + | |
427 | + | std::cerr << "ji[" << i << "]\t" |
428 | + | << ji[0] << "\t" |
429 | + | << ji[1] << "\t" |
430 | + | << ji[2] << "\n"; |
431 | } | |
432 | } | |
433 | ||
# | Line 690 | Line 683 | void Integrator::rotate( int axes1, int axes2, double | |
683 | ||
684 | ||
685 | void Integrator::rotate( int axes1, int axes2, double angle, double ji[3], | |
686 | < | double A[3][3] ){ |
686 | > | double A[9] ){ |
687 | ||
688 | int i,j,k; | |
689 | double sinAngle; | |
# | Line 702 | Line 695 | void Integrator::rotate( int axes1, int axes2, double | |
695 | double tempA[3][3]; | |
696 | double tempJ[3]; | |
697 | ||
698 | + | |
699 | // initialize the tempA | |
700 | ||
701 | for(i=0; i<3; i++){ | |
702 | for(j=0; j<3; j++){ | |
703 | < | tempA[j][i] = A[i][j]; |
703 | > | tempA[j][i] = A[3*i+j]; |
704 | } | |
705 | } | |
706 | ||
# | Line 763 | Line 757 | void Integrator::rotate( int axes1, int axes2, double | |
757 | ||
758 | for(i=0; i<3; i++){ | |
759 | for(j=0; j<3; j++){ | |
760 | < | A[j][i] = 0.0; |
760 | > | A[3*j+i] = 0.0; |
761 | for(k=0; k<3; k++){ | |
762 | < | A[j][i] += tempA[i][k] * rot[j][k]; |
762 | > | A[3*j+i] += tempA[i][k] * rot[j][k]; |
763 | } | |
764 | } | |
765 | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |