39c39
<   {ex0=NULL; mQQt=NULL; mQ=NULL; mhQ=NULL; /*tmpFX=NULL;*/ mhQhQt=NULL; tmpC=NULL; invSqrtMass=NULL;}
---
> 		{ex0=NULL; mQQt=NULL; mQ=NULL; mhQ=NULL; tmpFX=NULL; mhQhQt=NULL; tmpC=NULL; invSqrtMass=NULL;}
45c45
<   {ex0=NULL; mQQt=NULL; mQ=NULL; mhQ=NULL; /*tmpFX=NULL;*/ mhQhQt=NULL; tmpC=NULL; invSqrtMass=NULL;}
---
> 		{ex0=NULL; mQQt=NULL; mQ=NULL; mhQ=NULL; tmpFX=NULL; mhQhQt=NULL; tmpC=NULL; invSqrtMass=NULL;}
54c54
<     //    if(tmpFX!=NULL) delete [] tmpFX;
---
> 	  if(tmpFX!=NULL) delete [] tmpFX;
71c71
<     (*myPositions) -= cofm; 
---
> 	for( unsigned int i = 0; i < myPositions->size(); i++ ) (*myPositions)[i] -= cofm;
75c75
<     ex0->intoAssign(*myPositions);
---
>     *ex0 = *myPositions;
107,109c107
<       for (int j=0;j<_rfM;j++) mhQ[i + j*_3N] = eigvectp->myData[3*(i/3 + j*_N)+(i%3)];
<       //for(int j=0;j<_rfM;j++) mhQ[i + j*_3N] = (*eigvectp)[i/3 + j*_N][i%3];
<       //_3N
---
> 		for(int j=0;j<_rfM;j++) mhQ[i + j*_3N] = (*eigvectp)[i/3 + j*_N][i%3];//_3N
150c148
<     //    tmpFX = new double[_3N];
---
> 	tmpFX = new double[_3N];
167d164
< 
176a174
> 
180,182c178,181
<   //  double* NormModeInt::vector3DBlockTOvect(Vector3DBlock* blkDat){
<   //  return blkDat->myData;
<   //  }
---
>   double* NormModeInt::vector3DBlockTOvect(Vector3DBlock* blkDat, double* vecDat){
> 	  for( int i=0; i<3*_N; i++) vecDat[i] = (*blkDat)[i/3][i%3];
> 	  return vecDat;
>   }
185,188c184,187
<   //Vector3DBlock* NormModeInt::vectTOvector3DBlock(double* vecDat, Vector3DBlock* blkDat){
<   //  blkDat->myData = vecDat;
<   //  return blkDat;
<   // }
---
>   Vector3DBlock* NormModeInt::vectTOvector3DBlock(double* vecDat, Vector3DBlock* blkDat){
> 	  for( int i=0; i<3*_N; i++) (*blkDat)[i/3][i%3] = vecDat[i];
> 	  return blkDat;
>   }
193,194c192
< 	//tmpFX = vector3DBlockTOvect(force);	
< 	//put positions-x_0 into linear array
---
> 	vector3DBlockTOvect(force, tmpFX);	//put positions-x_0 into linear array
204c202
< 	dgemv_ (*transA, m, n, alpha, mhQ, m, force->myData, incxy, beta, tmpC, incxy, len_transa);
---
> 	dgemv_ (*transA, m, n, alpha, mhQ, m, tmpFX, incxy, beta, tmpC, incxy, len_transa);
212c210
< 	dgemv_ (*transB, m, n, alpha, mhQ, m, tmpC, incxy, beta, iPforce->myData, incxy, len_transa);
---
> 	dgemv_ (*transB, m, n, alpha, mhQ, m, tmpC, incxy, beta, tmpFX, incxy, len_transa);
216c214
< 	//vectTOvector3DBlock(tmpFX, iPforce);
---
> 	vectTOvector3DBlock(tmpFX, iPforce);
224,225c222
< 	//tmpFX = vector3DBlockTOvect(force);	
< 	//put positions-x_0 into linear array
---
> 	vector3DBlockTOvect(force, tmpFX);	//put positions-x_0 into linear array
235c232
< 	dgemv_ (*transA, m, n, alpha, mQ, m, force->myData, incxy, beta, tmpC, incxy, len_transa);
---
> 	dgemv_ (*transA, m, n, alpha, mQ, m, tmpFX, incxy, beta, tmpC, incxy, len_transa);
243c240
< 	dgemv_ (*transB, m, n, alpha, mQ, m, tmpC, incxy, beta, iPforce->myData, incxy, len_transa);
---
> 	dgemv_ (*transB, m, n, alpha, mQ, m, tmpC, incxy, beta, tmpFX, incxy, len_transa);
247c244
< 	//vectTOvector3DBlock(tmpFX, iPforce);
---
> 	vectTOvector3DBlock(tmpFX, iPforce);
255d251
<     Real invH = 1/h;
258c254
<     mySsVelocities.intoAssign(*myVelocities);
---
> 	mySsVelocities=*myVelocities;
270c266,269
<       *myVelocities += myForces->invMassScale(myTopo->atoms)*h*0.5;
---
> 	  //subspaceForce(myForces, myForces);	//moved!
> 	  //doHalfKick();
> 	  for( unsigned int i = 0; i < count; ++i)
> 			(*myVelocities)[i] += (*myForces)[i] * h * 0.5 / myTopo->atoms[i].scaledMass;
271a271
> 	  //if(!myNVE) subspaceForce(myVelocities, myVelocities);
275c275,277
< 	*myPositions += (*myVelocities) * h;
---
> 		//subspaceForce(myVelocities, &mySsVelocities);//myVelocities);//
> 		for( unsigned int i = 0; i < count; ++i)
> 			(*myPositions)[i]  += (*myVelocities)[i] * h;//mySsVelocities[i] * h;//
281c283
< 	oldPos.intoAssign(*myPositions);
---
> 		oldPos = *myPositions;
283c285,286
< 	  *myPositions += compVelocities*h;
---
> 			for( unsigned int i = 0; i < count; ++i)
> 				(*myPositions)[i] += compVelocities[i] * h;
290c293,296
< 	compVelocities = (*myPositions - oldPos)*invH;
---
> 		for( unsigned int i = 0; i < count; ++i){
> 			compVelocities[i] = ((*myPositions)[i] - oldPos[i]) / h;
> 			//(*myVelocities)[i] += ((*myPositions)[i] - oldPos[i] ) / h;
> 		}
293a300
> 	  //for( unsigned int i = 0; i < _N; ++i)	(*myForces)[i] /= myTopo->atoms[i].scaledMass;
295c302,307
<       *myVelocities += myForces->invMassScale(myTopo->atoms)*h*0.5;
---
> 	  //for( unsigned int i = 0; i < _N; ++i)	(*myForces)[i] *= myTopo->atoms[i].scaledMass;
> 	  //subspaceForce(myForces, myForces);
> 	  //
> 	  for( unsigned int i = 0; i < count; ++i)
> 			(*myVelocities)[i] += (*myForces)[i] * h * 0.5 / myTopo->atoms[i].scaledMass;
> 	  //doHalfKick();
296a309
> 	  //if(!myNVE) subspaceForce(myVelocities, myVelocities);
307c320
< 	*myVelocities *= efact;
---
> 		  for( int k = 0; k < _N; k++ ) (*myVelocities)[k] *= efact;
377,379c390,391
< 	unsigned int index = i*3+j;
< 	gaussRandCoord1[index] = randomGaussianNumber(mySeed) / sqrtmass;
< 	gaussRandCoord2[index] = randomGaussianNumber(mySeed) / sqrtmass;
---
> 			gaussRandCoord1[i][j] = randomGaussianNumber(mySeed) / sqrtmass;
> 			gaussRandCoord2[i][j] = randomGaussianNumber(mySeed) / sqrtmass;
385d396
<     //
388a400
> 	//
390c402
<     dgemv_ (*transA, m, n, alpha, mQ, m, tmpC, incxy, beta, gaussRandCoord1.myData, incxy, len_transa);
---
> 	dgemv_ (*transA, m, n, alpha, mQ, m, tmpC, incxy, beta, tmpFX, incxy, len_transa);
395c407
<       for( unsigned int j=0; j < 3; j++)  gaussRandCoord1.myData[3*i + j] *= invSqrtMass[i];//oosqrtmass;//
---
> 		for( unsigned int j=0; j < 3; j++)  tmpFX[3*i + j] *= invSqrtMass[i];//oosqrtmass;//
398c410
<     //vectTOvector3DBlock(tmpFX, &gaussRandCoord1);
---
> 	vectTOvector3DBlock(tmpFX, &gaussRandCoord1);
403c415
<     dgemv_ (*transA, m, n, alpha, mQ, m, tmpC, incxy, beta, gaussRandCoord2.myData, incxy, len_transa);
---
> 	dgemv_ (*transA, m, n, alpha, mQ, m, tmpC, incxy, beta, tmpFX, incxy, len_transa);
405,408d416
<     // Update first with Gauss1
<     posUpdt = (gaussRandCoord1*langDriftZ1 + (*myVelocities))*tau1;
<     *myVelocities = ((*myVelocities)*exp(-myGamma*dt))+(gaussRandCoord1*sqrtFCoverM*sqrtTau2);
< 
412c420
<       for( unsigned int j=0; j < 3; j++)  gaussRandCoord2.myData[3*i + j] *= invSqrtMass[i];//oosqrtmass;//
---
> 		for( unsigned int j=0; j < 3; j++)  tmpFX[3*i + j] *= invSqrtMass[i];//oosqrtmass;//
415c423
<     //vectTOvector3DBlock(tmpFX, &gaussRandCoord2);
---
> 	vectTOvector3DBlock(tmpFX, &gaussRandCoord2);
422,425c430,442
<     posUpdt += gaussRandCoord2*langDriftZ2*tau1;
<     //for( unsigned int i = 0; i < myPositions->size(); i++ ) {
<     //posUpdt[i] =(gaussRandCoord1[i]*langDriftZ1 +gaussRandCoord2[i]*langDriftZ2 +(*myVelocities)[i])*tau1;
<     
---
>     for( unsigned int i = 0; i < myPositions->size(); i++ ) {
> 	  //Real sqrtmass = sqrt( myTopo->atoms[i].scaledMass );
>       //  Generate gaussian random numbers for each spatial directions such that 
>       //  the average is 0 and the standard deviation is fParam.  The algorithm 
>       //  that is  used was taken from X-PLOR. The algorithm is a "sum of 
>       //  uniform deviates algorithm" which may be found in Abramowitz and 
>       //  Stegun, "Handbook of Mathematical Functions", pg 952.
>       // update drift(fluctuation)
> 	  //####FIX this!! wont be in subspace!
>       //(*myPositions)[i]+=(gaussRandCoord1[i]*langDriftZ1 / sqrtmass+gaussRandCoord2[i]*langDriftZ2 / sqrtmass+(*myVelocities)[i])*tau1;
> 	  //
> 	  posUpdt[i] =(gaussRandCoord1[i]*langDriftZ1 +gaussRandCoord2[i]*langDriftZ2 +(*myVelocities)[i])*tau1;
> 	  //posUpdt[i] =(gaussRandCoord1[i]*langDriftZ1 / sqrtmass+gaussRandCoord2[i]*langDriftZ2 / sqrtmass+(*myVelocities)[i])*tau1;
427,429c444,446
<     //(*myVelocities)[i] = (*myVelocities)[i]*exp(-myGamma*dt)+gaussRandCoord1[i]*sqrtFCoverM*sqrtTau2;
<     
<     //}
---
>       (*myVelocities)[i] = (*myVelocities)[i]*exp(-myGamma*dt)+gaussRandCoord1[i]*sqrtFCoverM*sqrtTau2;
>       //(*myVelocities)[i] = (*myVelocities)[i]*exp(-myGamma*dt)+gaussRandCoord1[i]*sqrtFCoverM*sqrtTau2 / sqrtmass;
>     }
432c449,450
<     *myPositions += posUpdt;
---
>     for( unsigned int i = 0; i < myPositions->size(); i++ )
> 			(*myPositions)[i]+=posUpdt[i];
435a454,455
> 
>