一个温控方法

 double cc,sumv[4],sumvc0,ct;
 int i,j,k;
 ek0=t*kb*node*1.5;
 cc=dt/mass;

     for(j=1;j<=3;j++)
  {
      sumv[j]=0.0;
   for(i=1;i<=node;i++)
   {
   
    temp0[i][j]=v[i][j]+f[i][j]*cc;
    sumv[j]+=temp0[i][j]*v[i][j];
   }     
  
  }
 for(i=1;i<=3;i++)
 {
  sumvc0+=sumv[i];
 }

 ct=(1-sqrt(sumvc0*mass/2./ek0))/dt;

 for(j=1;j<=3;j++)
 {
  for(i=1;i<=node;i++)
  {
  
    temp[i][j]=f[i][j]/mass+ct*v[i][j];
    v[i][j]=v[i][j]+dt*temp[i][j];
  }
 }

 for(j=1;j<=3;j++)
 {
  for(i=1;i<=node;i++)
  {
      
    r[i][j]=r[i][j]+dt*v[i][j];
  }
 }