70d69
<     Real* forcedata = myForces->myData;
78,79d76
< 	Real* pkdata = pk.myData;
< 	Real* gkdata = gk.myData;
82c79
< 	for(int i=0; i<_3N; i++) pkdata[i] = 0.0;
---
> 	for(int i=0; i<_3N; i++) pk[i/3][i%3] = 0.0;
87d83
< 	  //cout << "IN: " << in << " " << myEnergies->potentialEnergy() << endl;
95d90
< 	    cout << "ZEO" << endl;
98,102c93
< 	      if(lup>0) {
< 		cout << "CALCULATING FORCES..." << endl;
< 		calculateForces();
< 		cout << "DONE AGAIN" << endl;
< 	      }
---
> 				if(lup>0) calculateForces();
105,109c96
< 	      for( int k = 0; k < _3N; k += 3 )
< 		lambdaSlp -= pkdata[k]*forcedata[k] +
< 		  pkdata[k+1]*forcedata[k+1] +
< 		  pkdata[k+2]*forcedata[k+2];
< 	      
---
> 				for( int k = 0; k < _N; k++ ) lambdaSlp -= pk[k].dot((*myForces)[k]);
123c110
< 	      myPositions->intoAssign(oldPos);
---
> 				*myPositions = oldPos;
128d114
< 
136c122
< 	  oldGk.intoAssign(gk);
---
> 		oldGk = gk;
140a127
> 
142d128
< 	  Real* oldgkdat = oldGk.myData;
144,146c130,131
< 	    unsigned int index = k*3;
< 	    betaD += oldgkdat[index]*oldgkdat[index]+oldgkdat[index+1]*oldgkdat[index+1]+oldgkdat[index+2]*oldgkdat[index+2];
< 	    betaN += gkdata[index]*(gkdata[index]-oldgkdat[index])+gkdata[index+1]*(gkdata[index+1]-oldgkdat[index+1])+gkdata[index+2]*(gkdata[index+2]-oldgkdat[index+2]);
---
> 			betaD += oldGk[k].dot(oldGk[k]);//pk[k].dot(gk[k] - oldGk[k]);//
> 			betaN += gk[k].dot(gk[k] - oldGk[k]);
152d136
< 
156d139
< 
166,167c149
< 
< 	  //myPrevMTS->tmpFX = myPrevMTS->vector3DBlockTOvect(&pk);	
---
> 		myPrevMTS->vector3DBlockTOvect(&pk, myPrevMTS->tmpFX);	
170d151
< 
172,175c153
< 	  //dgemv_ (*transA, m, n, alpha, &myPrevMTS->mhQ[_3N*myPrevMTS->_rfM], m, myPrevMTS->tmpFX, incxy, beta, myPrevMTS->tmpC, incxy, len_transa);
< 
< 	  dgemv_ (*transA, m, n, alpha, &myPrevMTS->mhQ[_3N*myPrevMTS->_rfM], m, pkdata, incxy, beta, myPrevMTS->tmpC, incxy, len_transa);
< 
---
> 		dgemv_ (*transA, m, n, alpha, &myPrevMTS->mhQ[_3N*myPrevMTS->_rfM], m, myPrevMTS->tmpFX, incxy, beta, myPrevMTS->tmpC, incxy, len_transa);
182,183c160
< 	  //dgemv_ (*transB, m, n, alpha, &myPrevMTS->mhQ[_3N*myPrevMTS->_rfM], m, myPrevMTS->tmpC, incxy, beta, myPrevMTS->tmpFX, incxy, len_transa);
< 	  dgemv_ (*transB, m, n, alpha, &myPrevMTS->mhQ[_3N*myPrevMTS->_rfM], m, myPrevMTS->tmpC, incxy, beta, pkdata, incxy, len_transa);
---
> 		dgemv_ (*transB, m, n, alpha, &myPrevMTS->mhQ[_3N*myPrevMTS->_rfM], m, myPrevMTS->tmpC, incxy, beta, myPrevMTS->tmpFX, incxy, len_transa);
186c163
< 	  //myPrevMTS->vectTOvector3DBlock(myPrevMTS->tmpFX, &pk);
---
> 		myPrevMTS->vectTOvector3DBlock(myPrevMTS->tmpFX, &pk);
192c169
< 	  oldPos.intoAssign(*myPositions);
---
> 		oldPos = *myPositions;
198,199c175
< 	    int index = 3*(i/3+(_m-1)*_N)+(i%3);
< 	    if(maxEv < fabs((*myPrevMTS->eigvectp).myData[index])) {
---
> 			if(maxEv < fabs((*myPrevMTS->eigvectp)[i/3+(_m-1)*_N][i%3])){
201c177
< 	      maxEv = fabs((*myPrevMTS->eigvectp).myData[index]);
---
> 				maxEv = fabs((*myPrevMTS->eigvectp)[i/3+(_m-1)*_N][i%3]);
205,206c181
< 	  lambda = 0.5 * maxEv / (*myPrevMTS->eigvalp).myData[3*((_m-1)/3)+((_m-1)%3)];
< 
---
> 		lambda = 0.5 * maxEv / (*myPrevMTS->eigvalp)[(_m-1)/3][(_m-1)%3]; //Factor of 2? 
209a185
> 
211,216c187,188
< 	  //for( int k = 0; k < _N; k++ ) lambdaSlp2 -= pk[k].dot((*myForces)[k]);
< 	  for (int k = 0; k < _3N; k += 3) {
< 	    lambdaSlp2 -= (pkdata[k]*forcedata[k] +
< 			   pkdata[k+1]*forcedata[k+1] +
< 			   pkdata[k+2]*forcedata[k+2]);
< 	  }
---
> 		for( int k = 0; k < _N; k++ ) lambdaSlp2 -= pk[k].dot((*myForces)[k]);
> 		//report <<debug(7)<<"[NormModeMin::updtConstrainedModes] lambdaSlp0= "<<lambdaSlp<<" lambdSL2= "<<1e-4*lambda*lambdaSlp<<endl;
219d190
< 
221,222d191
< 
< 	    
228,233c197
< 	  for( int k = 0; k < _N; k++ ) {
< 	    unsigned int index = k*3;
< 	    lambdaSlp -= pkdata[index]*forcedata[index];
< 	    lambdaSlp -= pkdata[index+1]*forcedata[index+1];
< 	    lambdaSlp -= pkdata[index+2]*forcedata[index+2];
< 	  }
---
> 		for( int k = 0; k < _N; k++ ) lambdaSlp -= pk[k].dot((*myForces)[k]);
242c206
< 	  myPositions->intoAssign(oldPos);
---
> 		*myPositions = oldPos;
248,253c212
< 	  for( int k = 0; k < _N; k++ ) {
< 	    unsigned int index = k*3;
< 	    lambdaSlp -= pkdata[index]*forcedata[index];
< 	    lambdaSlp -= pkdata[index+1]*forcedata[index+1];
< 	    lambdaSlp -= pkdata[index+2]*forcedata[index+2];
< 	  }
---
> 		for( int k = 0; k < _N; k++ ) lambdaSlp -= pk[k].dot((*myForces)[k]);
287,293c246
<     Real* thedata = blk2norm.myData;
<     for( int k = 0; k < _N; k++ ) {
<       int index = k*3;
<       normTot += (thedata[index]*thedata[index] +
< 		  thedata[index+1]*thedata[index+1] +
< 		  thedata[index+2]*thedata[index+2]);
<     }
---
> 	for( int k = 0; k < _N; k++ ) normTot += blk2norm[k].normSquared();// / myTopo->atoms[k].scaledMass;