Nose-Hoover控温

 double sumv,nita,tao;
 int i,j,k;

 sumv=0.0;
 for(i=1;i<=node;i++)
 {
  for(j=1;j<=3;j++)
  {
  
    r[i][j]=r[i][j]+dt*v[i][j]+0.5*dt*dt*(f[i][j]/mass-nita*v[i][j]);
    sumv+=v[i][j]*v[i][j];
  }
 }
 t0=mass*sumv/(3*node*kb);
    nita=nita+0.5*dt*(t0-t)/t/(tao*tao);
 sumv=0.0;
 for(i=1;i<=node;i++)
 {
  for(j=1;j<=3;j++)
  {
  
    v[i][j]=v[i][j]+0.5*dt*(f[i][j]/mass-nita*v[i][j]);
    sumv+=v[i][j]*v[i][j]
  }
 }
    t0=mass*sumv/(3*node*kb);
    nita=nita+0.5*dt*(t0-t)/t/(tao*tao);

 list();
 force();
 
 for(i=1;i<=node;i++)
 {
  for(j=1;j<=3;j++)
  {
  
    v[i][j]=1/(1+0.5*dt*nita)*(v[i][j]+0.5*dt*f[i][j]/mass);

  }
 }