Permission denied, please try again.
Index: LangevinImpulseIntegrator.cpp
===================================================================
RCS file: /cvsroot/protomol/protomol/framework/integrators/LangevinImpulseIntegrator.cpp,v
retrieving revision 1.13
diff -r1.13 LangevinImpulseIntegrator.cpp
42a43
>     cout << "STSS INIT" << endl;
43a45
>     cout << "FORCE INIT" << endl;
44a47
>     cout << "LALALA" << endl;
59a63
>     const Real tau3 = myGamma*dt;
71a76,97
> 
>        const Vector3DBlock *forces = getForces();
>       
>        Real wp, wn;
>        wp = (tau1-1+tau3)/(tau3*(1-tau1));
>        if(thisStep >= (totalSteps-1)) wn = 1-wp;
>        else wn = 0;
>      
> 
>      Vector3DBlock grn1;
> 
>      bool betaFlag;
> 
>      if(beta.size() != myTopo->atoms.size()) betaFlag = true;
>      else betaFlag = false; 
> 
>      Real* veldata = myVelocities->myData;
>      Real* posdata = myPositions->myData;
>      Real* forcedata = myForces->myData;
> 
> //     for( unsigned int i = 0; i < myVelocities->size(); i++ ) {
> // =======
74,85d99
<       Real sqrtFCoverM = sqrt( forceConstant / mass );
<       Real langDriftVal = sqrtFCoverM / myGamma;
<       Real langDriftZ1 = langDriftVal * ( tau1 - tau2 ) / sqrtTau2;
<       Real langDriftZ2 = langDriftVal * sqrtVal1;
< 
<       //  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.
<       Vector3D gaussRandCoord1(randomGaussianNumber(mySeed),randomGaussianNumber(mySeed),randomGaussianNumber(mySeed));
<       Vector3D gaussRandCoord2(randomGaussianNumber(mySeed),randomGaussianNumber(mySeed),randomGaussianNumber(mySeed));
87,88c101,163
<       // update drift(fluctuation)
<       (*myPositions)[i]+=(gaussRandCoord1*langDriftZ1+gaussRandCoord2*langDriftZ2+(*myVelocities)[i])*tau1;
---
> 
>       Real ia,ib,ic;
>       ia = a/mass;
>       ib = b/mass;
>       ic = c/mass;
> 
> 
>       if(thisStep == totalSteps) {
>           beta[i] = ib/alpha[i];
>           alpha[i] = sqrt(ic-pow(beta[i],2));
>       }
>       else {
> 	  if(betaFlag) {
>           	beta.push_back(ib/alpha[i]);
>           	alpha[i] = sqrt(ia+ic-pow(beta[i],2));
> 	  }else {
>         	beta[i] = ib/alpha[i];
>         	alpha[i] = sqrt(ia+ic-pow(beta[i],2));
>       	}
>       }
>         
>       Vector3D gaussRandCoord(randomGaussianNumber(mySeed),randomGaussianNumber(mySeed),randomGaussianNumber(mySeed));
> 
>       Vector3D Rn(gaussRandn1.myData[i*3]*beta[i]+gaussRandCoord.x*alpha[i], gaussRandn1.myData[i*3+1]*beta[i]+gaussRandCoord.y*alpha[i], gaussRandn1.myData[i*3+2]*beta[i]+gaussRandCoord.z*alpha[i]);
> 
>       int index = i*3;
> 
>       if(thisStep == totalSteps) {
> 	//full update of velocities
> 	veldata[index] = (veldata[index])*tau2+(forcedata[index])*dt*wn*(1/mass)+Rn.x;
> 	veldata[index+1] = (veldata[index+1])*tau2+(forcedata[index+1])*dt*wn*(1/mass)+Rn.y;
> 	veldata[index+2] = (veldata[index+2])*tau2+(forcedata[index+2])*dt*wn*(1/mass)+Rn.z;
> 	continue;
>       }
> 
>       
>       //update velocity at half step
>       veldata[index] += (veldata[index]*tau2+(forcedata[index])*dt*(1/mass)+Rn.x)*tau2;
>       veldata[index+1] += (veldata[index+1]*tau2+(forcedata[index+1])*dt*(1/mass)+Rn.y)*tau2;
>       veldata[index+2] += (veldata[index+2]*tau2+(forcedata[index+2])*dt*(1/mass)+Rn.z)*tau2;
> 
>       //update positions at full step
>       posdata[index] += veldata[index]*((1-tau1)/(myGamma*tau2));
>       posdata[index+1] += veldata[index+1]*((1-tau1)/(myGamma*tau2));
>       posdata[index+2] += veldata[index+2]*((1-tau1)/(myGamma*tau2));
> 
>       grn1.push_back(gaussRandCoord);
> // =======
> //       Real sqrtFCoverM = sqrt( forceConstant / mass );
> //       Real langDriftVal = sqrtFCoverM / myGamma;
> //       Real langDriftZ1 = langDriftVal * ( tau1 - tau2 ) / sqrtTau2;
> //       Real langDriftZ2 = langDriftVal * sqrtVal1;
> 
> //       //  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.
> //       Vector3D gaussRandCoord1(randomGaussianNumber(mySeed),randomGaussianNumber(mySeed),randomGaussianNumber(mySeed));
> //       Vector3D gaussRandCoord2(randomGaussianNumber(mySeed),randomGaussianNumber(mySeed),randomGaussianNumber(mySeed));
> 
> //       // update drift(fluctuation)
> //       (*myPositions)[i]+=(gaussRandCoord1*langDriftZ1+gaussRandCoord2*langDriftZ2+(*myVelocities)[i])*tau1;
90,91c165,167
<       // semi-update velocities
<       (*myVelocities)[i] = (*myVelocities)[i]*exp(-myGamma*dt)+gaussRandCoord1*sqrtFCoverM*sqrtTau2;
---
> //       // semi-update velocities
> //       (*myVelocities)[i] = (*myVelocities)[i]*exp(-myGamma*dt)+gaussRandCoord1*sqrtFCoverM*sqrtTau2;
> // >>>>>>> 1.13
93a170,220
> 
> //   void LangevinImpulseIntegrator::startFluctuationKick() {
> 
> //     Real* veldata = myVelocities->myData;
> //     Real* posdata = myPositions->myData;
> //     Real* forcedata = myForces->myData;
> 
> 
> //       const Vector3DBlock *forces = getForces(); 
> //       const Real dt = getTimestep() * Constant::INV_TIMEFACTOR; 
> //       myGamma = myGamma/Constant::INV_TIMEFACTOR;
> //       const Real tau1 = exp( -myGamma*dt);
> //       const Real tau2 = exp( -myGamma*(dt/2));
> //       const Real tau3 = myGamma*dt;
> 
> //       Real wp = (tau1-1+tau3)/(tau3*(1-tau1));
> //       Real wn = 1 -wp;
> 
> //        a = Constant::BOLTZMANN*myLangevinTemperature*(2*pow(wp,2)*myGamma*dt+wp-wn);
> //        b = Constant::BOLTZMANN*myLangevinTemperature*(2*wp*wn*myGamma*dt+wn-wp);
> //        c = Constant::BOLTZMANN*myLangevinTemperature*(2*pow(wn,2)*myGamma*dt+wp-wn);
> 
> 
> //      Vector3DBlock grn1; 
> 
> //     for( unsigned int i = 0; i < myVelocities->size(); i++ ) {
> //       int index = i*3;
> //       Real mass = myTopo->atoms[i].scaledMass;
> //       Real ia;
> //       ia = a/mass;
> //       alpha.push_back(sqrt(ia));
> //       Vector3D gaussRandCoord(randomGaussianNumber(mySeed),randomGaussianNumber(mySeed),randomGaussianNumber(mySeed));
> //       Vector3D R0 = gaussRandCoord*alpha[i];
> 
> //       //Update of half step velocities
> //       veldata[index] += (forcedata[index]*dt*wp*(1/mass)+R0.x)*tau2;
> //       veldata[index+1] += (forcedata[index+1]*dt*wp*(1/mass)+R0.y)*tau2;
> //       veldata[index+2] += (forcedata[index+2]*dt*wp*(1/mass)+R0.z)*tau2;
> 
> //       //Updating full step positions
> //       posdata[index] += veldata[index]*((1-tau1)/(myGamma*tau2));
> //       posdata[index+1] += veldata[index+1]*((1-tau1)/(myGamma*tau2));
> //       posdata[index+2] += veldata[index+2]*((1-tau1)/(myGamma*tau2));
> 
> //       grn1.push_back(gaussRandCoord);
> //     }
> //     gaussRandn1 = grn1;
>       
> //   }  
> 
>