Calculation of Energy Spectrum using C++
Avoid using namespace std;
If you are coding professionally you probably should get out of the habit of using the using namespace std; statement. The code will more clearly define where cout and other identifiers are coming from (std::cin, std::cout). As you start using namespaces in your code it is better to identify where each function comes from because there may be function name collisions from different namespaces. The identifiercout you may override within your own classes, and you may override the operator &lt;&lt; in your own classes as well. This stack overflow question discusses this in more detail.
Include Files
You are currently including math.h 3 times, it is only necessary to do this once. Since you are using C++, stick with #include &lt;cmath&gt;.
Declare the Variables as Needed
The code contains this line:
    int i,j,k;

These are the variables you are using as loop control variables, they can therefore be declared in the loops themselves:
    for (int i = 0; i &lt; np-1; i++){
        for (int j = i+1; j &lt; np; j++){
        ...
        }
    }

Break the Code Into Functions to Make it Easier to Implement and Debug
Currently the code is all in main(), as programs grow and features get added the code would be to complex to easily update if it is all in main(). Also to achieve the ultimate goal of multi-threading for this program you will need to break the code up into multiple functions that can run on a separate core.
void FUNCTION_NAME(int i, int j, int np, double E[])
{
    double x = sin(double(2*i+j)/(3.0*np));//x,y,z,d1 for test.
    double y = sin(double(i+j)/(2.0*np));//x,y,z are velocity
    double z = sin(double(i*j)/np/np); //d1 velocity magnitude
    double d1 = sqrt(x*x + y*y + z*z);

    for (k = 1; k &lt;= kmax; k++){
        double d2 = k * d1;
        E[k] += x + y + sin(d2) / d2;
    }
}

int main()
{
    int np = 1000000;  // number of fluid particle
    int kmax = 100;   // wavenumber
    double E[kmax+1] = {0.0};  //energy

    // need to be parallelized
    for (i = 0; i &lt; np-1; i++){
        for (j = i+1; j &lt; np; j++){
            FUNCTION_NAME(i, j, np, E[]);
        }
    }
}